SphericalVertexCluster.hpp

Go to the documentation of this file.

#pragma once
 
#include <Eigen/Core>
 
#include <boost/container/flat_set.hpp>
 
#include "mapping/RadialBasisFctSolver.hpp"
#include "mapping/config/MappingConfiguration.hpp"
#include "mapping/impl/BasisFunctions.hpp"
#include "mesh/Filter.hpp"
#include "precice/impl/Types.hpp"
#include "profiling/Event.hpp"
 
namespace precice::mapping {

template <typename RADIAL_BASIS_FUNCTION_T>
class SphericalVertexCluster {
public:
  SphericalVertexCluster(mesh::Vertex            center,
                         double                  radius,
                         RADIAL_BASIS_FUNCTION_T function,
                         Polynomial              polynomial,
                         mesh::PtrMesh           inputMesh,
                         mesh::PtrMesh           outputMesh);

  void mapConservative(const time::Sample &inData, Eigen::VectorXd &outData) const;

  void computeCacheData(const Eigen::Ref<const Eigen::MatrixXd> &globalIn, Eigen::MatrixXd &polyOut, Eigen::MatrixXd &coefficientsOut) const;

  void mapConsistent(const time::Sample &inData, Eigen::VectorXd &outData) const;

  void setNormalizedWeight(double normalizedWeight, VertexID vertexID);
 
  void evaluateConservativeCache(Eigen::MatrixXd &epsilon, const Eigen::MatrixXd &Au, Eigen::Ref<Eigen::MatrixXd> out);

  double computeWeight(const mesh::Vertex &v) const;
 
  Eigen::VectorXd interpolateAt(const mesh::Vertex &v, const Eigen::MatrixXd &poly, const Eigen::MatrixXd &coeffs, const mesh::Mesh &inMesh) const;
  unsigned int getNumberOfInputVertices() const;

  std::array<double, 3> getCenterCoords() const;

  void clear();

  bool empty() const;
 
  void addWriteDataToCache(const mesh::Vertex &v, const Eigen::VectorXd &load, Eigen::MatrixXd &epsilon, Eigen::MatrixXd &Au,
                           const mesh::Mesh &inMesh);
 
  void initializeCacheData(Eigen::MatrixXd &polynomial, Eigen::MatrixXd &coeffs, const int nComponents);
 
private:
  mutable precice::logging::Logger _log{"mapping::SphericalVertexCluster"};

  mesh::Vertex _center;

  const double _radius;

  RadialBasisFctSolver<RADIAL_BASIS_FUNCTION_T> _rbfSolver;
 
  // Stores the global IDs of the vertices so that we can apply a binary
  // search in order to query specific objects. Here we have logarithmic
  // complexity for setting the weights (only performed once), and constant
  // complexity when traversing through the IDs (performed in each iteration).

  boost::container::flat_set<VertexID> _inputIDs;

  boost::container::flat_set<VertexID> _outputIDs;

  Eigen::VectorXd _normalizedWeights;

  Polynomial _polynomial;
 
  RADIAL_BASIS_FUNCTION_T _function;

  CompactPolynomialC2 _weightingFunction;

  bool _hasComputedMapping = false;
};
 
// --------------------------------------------------- HEADER IMPLEMENTATIONS
 
template <typename RADIAL_BASIS_FUNCTION_T>
SphericalVertexCluster<RADIAL_BASIS_FUNCTION_T>::SphericalVertexCluster(
    mesh::Vertex            center,
    double                  radius,
    RADIAL_BASIS_FUNCTION_T function,
    Polynomial              polynomial,
    mesh::PtrMesh           inputMesh,
    mesh::PtrMesh           outputMesh)
    : _center(center), _radius(radius), _polynomial(polynomial), _function(function), _weightingFunction(radius)
{
  PRECICE_TRACE(_center.getCoords(), _radius);
  precice::profiling::Event eq("map.pou.computeMapping.queryVertices");
  // Disable integrated polynomial, as it might cause locally singular matrices
  PRECICE_ASSERT(_polynomial != Polynomial::ON, "Integrated polynomial is not supported for partition of unity data mappings.");
 
  // Get vertices to be mapped
  // Subtract a safety margin to exclude the vertices at the edge
  auto outIDs = outputMesh->index().getVerticesInsideBox(center, radius - math::NUMERICAL_ZERO_DIFFERENCE);
  // Constructing the partition when we don't have evaluation points is pointless
  auto inIDs = inputMesh->index().getVerticesInsideBox(center, radius);
 
  // Transform the vector to the appropriate boost data structure
  // The IDs are sorted in the boost flat_set, hence, the function here has N log(N) complexity
  _inputIDs.insert(inIDs.begin(), inIDs.end());
  _outputIDs.insert(outIDs.begin(), outIDs.end());
  eq.stop();
  // If the cluster is empty, we return immediately
  if (empty()) {
    return;
  }
 
  PRECICE_DEBUG("SphericalVertexCluster input size: {}", inIDs.size());
  PRECICE_DEBUG("SphericalVertexCluster output size: {}", outIDs.size());
 
  // The polynomial system is underdetermined if inIDs.size() < dimension + 1. However, the dynamic adoption of the axis in the RBF solver
  // disables axis, if necessary. Hence, we don't disable the complete polynomial here for underdetermined systems. The case should anyway
  // only occur for an almost unreasonable small vertices-per-cluster configuration.
 
  // Construct the solver. Here, the constructor of the RadialBasisFctSolver computes already the decompositions etc, such that we can mark the
  // mapping in this cluster as computed (mostly for debugging purpose)
  std::vector<bool>         deadAxis(inputMesh->getDimensions(), false);
  precice::profiling::Event e("map.pou.computeMapping.rbfSolver");
  _rbfSolver          = RadialBasisFctSolver<RADIAL_BASIS_FUNCTION_T>{function, *inputMesh.get(), _inputIDs, *outputMesh.get(), _outputIDs, deadAxis, _polynomial};
  _hasComputedMapping = true;
}
 
template <typename RADIAL_BASIS_FUNCTION_T>
void SphericalVertexCluster<RADIAL_BASIS_FUNCTION_T>::setNormalizedWeight(double normalizedWeight, VertexID id)
{
  PRECICE_ASSERT(_outputIDs.size() > 0);
  PRECICE_ASSERT(_outputIDs.contains(id), id);
  PRECICE_ASSERT(normalizedWeight > 0);
 
  if (_normalizedWeights.size() == 0)
    _normalizedWeights.resize(_outputIDs.size());
 
  // The find method of boost flat_set comes with O(log(N)) complexity (the more expensive part here)
  auto localID = _outputIDs.index_of(_outputIDs.find(id));
 
  PRECICE_ASSERT(static_cast<Eigen::Index>(localID) < _normalizedWeights.size(), localID, _normalizedWeights.size());
  _normalizedWeights[localID] = normalizedWeight;
}
 
template <typename RADIAL_BASIS_FUNCTION_T>
void SphericalVertexCluster<RADIAL_BASIS_FUNCTION_T>::mapConservative(const time::Sample &inData, Eigen::VectorXd &outData) const
{
  // First, a few sanity checks. Empty partitions shouldn't be stored at all
  PRECICE_ASSERT(!empty());
  PRECICE_ASSERT(_hasComputedMapping);
  PRECICE_ASSERT(_normalizedWeights.size() == static_cast<Eigen::Index>(_outputIDs.size()));
 
  // Define an alias for data dimension in order to avoid ambiguity
  const unsigned int nComponents = inData.dataDims;
  const auto        &localInData = inData.values;
 
  // TODO: We can probably reduce the temporary allocations here
  Eigen::VectorXd in(_rbfSolver.getOutputSize());
 
  // Now we perform the data mapping component-wise
  for (unsigned int c = 0; c < nComponents; ++c) {
    // Step 1: extract the relevant input data from the global input data and store
    // it in a contiguous array, which is required for the RBF solver
    for (unsigned int i = 0; i < _outputIDs.size(); ++i) {
      const auto dataIndex = *(_outputIDs.nth(i));
      PRECICE_ASSERT(dataIndex * nComponents + c < localInData.size(), dataIndex * nComponents + c, localInData.size());
      PRECICE_ASSERT(_normalizedWeights[i] > 0, _normalizedWeights[i], i);
      // here, we also directly apply the weighting, i.e., we split the input data
      in[i] = localInData[dataIndex * nComponents + c] * _normalizedWeights[i];
    }
 
    // Step 2: solve the system using a conservative constraint
    auto result = _rbfSolver.solveConservative(in, _polynomial);
    PRECICE_ASSERT(result.size() == static_cast<Eigen::Index>(_inputIDs.size()));
 
    // Step 3: now accumulate the result into our global output data
    for (unsigned int i = 0; i < _inputIDs.size(); ++i) {
      const auto dataIndex = *(_inputIDs.nth(i));
      PRECICE_ASSERT(dataIndex * nComponents + c < outData.size(), dataIndex * nComponents + c, outData.size());
      outData[dataIndex * nComponents + c] += result(i);
    }
  }
}
 
template <typename RADIAL_BASIS_FUNCTION_T>
void SphericalVertexCluster<RADIAL_BASIS_FUNCTION_T>::evaluateConservativeCache(Eigen::MatrixXd &epsilon, const Eigen::MatrixXd &Au, Eigen::Ref<Eigen::MatrixXd> out)
{
  Eigen::MatrixXd localIn(_inputIDs.size(), Au.cols());
  _rbfSolver.evaluateConservativeCache(epsilon, Au, localIn);
  // Step 3: now accumulate the result into our global output data
  for (std::size_t i = 0; i < _inputIDs.size(); ++i) {
    const auto dataIndex = *(_inputIDs.nth(i));
    PRECICE_ASSERT(dataIndex < out.cols(), out.cols());
    out.col(dataIndex) += localIn.row(i);
  }
}
 
template <typename RADIAL_BASIS_FUNCTION_T>
void SphericalVertexCluster<RADIAL_BASIS_FUNCTION_T>::computeCacheData(const Eigen::Ref<const Eigen::MatrixXd> &globalIn, Eigen::MatrixXd &polyOut, Eigen::MatrixXd &coeffOut) const
{
  PRECICE_TRACE();
  Eigen::MatrixXd in(_rbfSolver.getInputSize(), globalIn.rows());
  // Step 1: extract the relevant input data from the global input data and store
  // it in a contiguous array, which is required for the RBF solver (last polyparams entries remain zero)
  for (std::size_t i = 0; i < _inputIDs.size(); i++) {
    const auto dataIndex = *(_inputIDs.nth(i));
    PRECICE_ASSERT(dataIndex < globalIn.cols(), globalIn.cols());
    in.row(i) = globalIn.col(dataIndex);
  }
  // Step 2: solve the system using a consistent constraint
  _rbfSolver.computeCacheData(in, _polynomial, polyOut, coeffOut);
}
 
template <typename RADIAL_BASIS_FUNCTION_T>
Eigen::VectorXd SphericalVertexCluster<RADIAL_BASIS_FUNCTION_T>::interpolateAt(const mesh::Vertex &v, const Eigen::MatrixXd &poly, const Eigen::MatrixXd &coeffs, const mesh::Mesh &inMesh) const
{
  PRECICE_TRACE();
  return _rbfSolver.interpolateAt(v, poly, coeffs, _function, _inputIDs, inMesh);
}
 
template <typename RADIAL_BASIS_FUNCTION_T>
void SphericalVertexCluster<RADIAL_BASIS_FUNCTION_T>::addWriteDataToCache(const mesh::Vertex &v, const Eigen::VectorXd &load,
                                                                          Eigen::MatrixXd &epsilon, Eigen::MatrixXd &Au,
                                                                          const mesh::Mesh &inMesh)
{
  PRECICE_TRACE();
  _rbfSolver.addWriteDataToCache(v, load, epsilon, Au, _function, _inputIDs, inMesh);
}
 
template <typename RADIAL_BASIS_FUNCTION_T>
void SphericalVertexCluster<RADIAL_BASIS_FUNCTION_T>::initializeCacheData(Eigen::MatrixXd &poly, Eigen::MatrixXd &coeffs, const int nComponents)
{
  if (Polynomial::SEPARATE == _polynomial) {
    poly.resize(_rbfSolver.getNumberOfPolynomials(), nComponents);
  }
  coeffs.resize(_inputIDs.size(), nComponents);
}
 
template <typename RADIAL_BASIS_FUNCTION_T>
void SphericalVertexCluster<RADIAL_BASIS_FUNCTION_T>::mapConsistent(const time::Sample &inData, Eigen::VectorXd &outData) const
{
  // First, a few sanity checks. Empty partitions shouldn't be stored at all
  PRECICE_ASSERT(!empty());
  PRECICE_ASSERT(_hasComputedMapping);
  PRECICE_ASSERT(_normalizedWeights.size() == static_cast<Eigen::Index>(_outputIDs.size()));
 
  // Define an alias for data dimension in order to avoid ambiguity
  const unsigned int nComponents = inData.dataDims;
  const auto        &localInData = inData.values;
 
  Eigen::VectorXd in(_rbfSolver.getInputSize());
 
  // Now we perform the data mapping component-wise
  for (unsigned int c = 0; c < nComponents; ++c) {
    // Step 1: extract the relevant input data from the global input data and store
    // it in a contiguous array, which is required for the RBF solver (last polyparams entries remain zero)
    for (unsigned int i = 0; i < _inputIDs.size(); i++) {
      const auto dataIndex = *(_inputIDs.nth(i));
      PRECICE_ASSERT(dataIndex * nComponents + c < localInData.size(), dataIndex * nComponents + c, localInData.size());
      in[i] = localInData[dataIndex * nComponents + c];
    }
 
    // Step 2: solve the system using a consistent constraint
    auto result = _rbfSolver.solveConsistent(in, _polynomial);
    PRECICE_ASSERT(static_cast<Eigen::Index>(_outputIDs.size()) == result.size());
 
    // Step 3: now accumulate the result into our global output data
    for (unsigned int i = 0; i < _outputIDs.size(); ++i) {
      const auto dataIndex = *(_outputIDs.nth(i));
      PRECICE_ASSERT(dataIndex * nComponents + c < outData.size(), dataIndex * nComponents + c, outData.size());
      PRECICE_ASSERT(_normalizedWeights[i] > 0);
      // here, we also directly apply the weighting, i.e., split the result data
      outData[dataIndex * nComponents + c] += result(i) * _normalizedWeights[i];
    }
  }
}
 
template <typename RADIAL_BASIS_FUNCTION_T>
double SphericalVertexCluster<RADIAL_BASIS_FUNCTION_T>::computeWeight(const mesh::Vertex &v) const
{
  // Assume that the local interpolant and the weighting function are the same
  // TODO: We don't need to reduce the dead coordinates here as the values should reduce anyway
  auto res = computeSquaredDifference(_center.rawCoords(), v.rawCoords(), {{true, true, true}});
  return _weightingFunction.evaluate(std::sqrt(res));
}
 
template <typename RADIAL_BASIS_FUNCTION_T>
unsigned int SphericalVertexCluster<RADIAL_BASIS_FUNCTION_T>::getNumberOfInputVertices() const
{
  return _inputIDs.size();
}
 
template <typename RADIAL_BASIS_FUNCTION_T>
std::array<double, 3> SphericalVertexCluster<RADIAL_BASIS_FUNCTION_T>::getCenterCoords() const
{
  return _center.rawCoords();
}
 
template <typename RADIAL_BASIS_FUNCTION_T>
bool SphericalVertexCluster<RADIAL_BASIS_FUNCTION_T>::empty() const
{
  return _inputIDs.size() == 0;
}
 
template <typename RADIAL_BASIS_FUNCTION_T>
void SphericalVertexCluster<RADIAL_BASIS_FUNCTION_T>::clear()
{
  PRECICE_TRACE();
  _inputIDs.clear();
  _outputIDs.clear();
  _rbfSolver.clear();
  _hasComputedMapping = false;
}
} // namespace precice::mapping