ITHACAutilities.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
     ██╗████████╗██╗  ██╗ █████╗  ██████╗ █████╗       ███████╗██╗   ██╗
     ██║╚══██╔══╝██║  ██║██╔══██╗██╔════╝██╔══██╗      ██╔════╝██║   ██║
     ██║   ██║   ███████║███████║██║     ███████║█████╗█████╗  ██║   ██║
     ██║   ██║   ██╔══██║██╔══██║██║     ██╔══██║╚════╝██╔══╝  ╚██╗ ██╔╝
     ██║   ██║   ██║  ██║██║  ██║╚██████╗██║  ██║      ██║      ╚████╔╝
     ╚═╝   ╚═╝   ╚═╝  ╚═╝╚═╝  ╚═╝ ╚═════╝╚═╝  ╚═╝      ╚═╝       ╚═══╝
 
 * In real Time Highly Advanced Computational Applications for Finite Volumes
 * Copyright (C) 2017 by the ITHACA-FV authors
-------------------------------------------------------------------------------
 
License
    This file is part of ITHACA-FV
 
    ITHACA-FV is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
 
    ITHACA-FV 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 Lesser General Public License for more details.
 
    You should have received a copy of the GNU Lesser General Public License
    along with ITHACA-FV. If not, see <http://www.gnu.org/licenses/>.
 
\*---------------------------------------------------------------------------*/
 
#include "ITHACAutilities.H"
#include "ITHACAstream.H"
#include "ITHACAparameters.H"
#include "turbulentTransportModel.H"
 
 
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
namespace ITHACAutilities
{
 
Eigen::MatrixXd rand(label rows, label cols, double min,
                     double max)
{
    std::srand(static_cast<long unsigned int>
               (std::chrono::high_resolution_clock::now().time_since_epoch().count()));
    Eigen::MatrixXd matr = Eigen::MatrixXd::Random(rows, cols);
    matr = (matr.array() + 1) / 2;
    matr = matr.array() * (max - min);
    matr = matr.array() + min;
    return matr;
}
 
Eigen::MatrixXd rand(label rows, Eigen::MatrixXd minMax)
{
    std::srand(static_cast<long unsigned int>
               (std::chrono::high_resolution_clock::now().time_since_epoch().count()));
    label cols = minMax.rows();
    Eigen::MatrixXd matr = Eigen::MatrixXd::Random(rows, cols);
    matr = (matr.array() + 1) / 2;
 
    for (label i = 0; i < cols; i++)
    {
        matr.col(i) = matr.col(i).array() * (minMax(i, 1) - minMax(i, 0));
        matr.col(i) = matr.col(i).array() + (minMax(i, 0));
    }
 
    return matr;
}
 
 
bool isInteger(double ratio)
{
    bool checkResult = 0;
 
    if (abs(round(ratio) - ratio) < std::sqrt(SMALL))
    {
        checkResult = true;
    }
    else
    {
        checkResult = false;
    }
 
    return checkResult;
}
 
bool isTurbulent()
{
    bool checkTurb;
    ITHACAparameters* para = ITHACAparameters::getInstance();
    auto& tur =
        para->mesh.lookupObject<incompressible::turbulenceModel>("turbulenceProperties");
 
    if (tur.type() == "Stokes" || tur.type() == "Maxwell"
            || tur.type() == "laminarModel")
    {
        checkTurb = false;
    }
    else
    {
        checkTurb = true;
    }
 
    return checkTurb;
}
 
template<typename T>
List<T> combineList(List<List<T>>& doubleList)
{
    List<T> a = ListListOps::combine<List<T>>(doubleList,
                accessOp<List<T>>());
#if OPENFOAM >= 1812
    inplaceUniqueSort(a);
#else
    labelList order;
    uniqueOrder(a, order);
    List<T> b(order.size());
 
    for (label i = 0; i < order.size(); ++i)
    {
        b[i] = a[order[i]];
    }
 
    a.resize(order.size());
    a = b;
#endif
    return a;
}
 
template List<label> combineList(List<List<label>>& doubleList);
 
 
// Using the Eigen library, using the SVD decomposition method to solve the
// matrix pseudo-inverse, the default error er is 0
Eigen::MatrixXd pinv_eigen_based(Eigen::MatrixXd &origin, const float er) {
  // perform svd decomposition
  Eigen::JacobiSVD<Eigen::MatrixXd> svd_holder(origin, Eigen::ComputeThinU | Eigen::ComputeThinV);
  // Build SVD decomposition results
  Eigen::MatrixXd U = svd_holder.matrixU();
  Eigen::MatrixXd V = svd_holder.matrixV();
  Eigen::MatrixXd D = svd_holder.singularValues();
 
  // Build the S matrix
  Eigen::MatrixXd S(V.cols(), U.cols());
  S.setZero();
 
  for (unsigned int i = 0; i < D.size(); ++i) {
 
    if (D(i, 0) > er) {
      S(i, i) = 1 / D(i, 0);
    } else {
      S(i, i) = 0;
    }
  }
  return V * S * U.transpose();
}
 
 
Eigen::MatrixXd invertMatrix(Eigen::MatrixXd& matrixToInvert, const word inversionMethod){
 
  Info << "Inversion method : " << inversionMethod << endl;
  if(inversionMethod == "pinv_eigen_based"){
    return pinv_eigen_based(matrixToInvert);
  }
  else if(inversionMethod == "direct"){
    return matrixToInvert.inverse();
  }
  else if(inversionMethod == "fullPivLu"){
    return matrixToInvert.fullPivLu().inverse();
  }
  else if(inversionMethod == "partialPivLu"){
    return matrixToInvert.partialPivLu().inverse();
  }
  else if(inversionMethod == "householderQr"){
    return matrixToInvert.householderQr().solve(Eigen::MatrixXd::Identity(matrixToInvert.rows(), matrixToInvert.cols()));
  }
  else if(inversionMethod == "colPivHouseholderQr"){
    return matrixToInvert.colPivHouseholderQr().inverse();
  }
  else if(inversionMethod == "fullPivHouseholderQr"){
    return matrixToInvert.fullPivHouseholderQr().inverse();
  }
  else if(inversionMethod == "completeOrthogonalDecomposition"){
    return matrixToInvert.completeOrthogonalDecomposition().pseudoInverse();
  }
  else if(inversionMethod == "jacobiSvd")
  {
    return matrixToInvert.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(Eigen::MatrixXd::Identity(matrixToInvert.rows(), matrixToInvert.cols()));
  }
  else if(inversionMethod == "llt")
  {
    return matrixToInvert.llt().solve(Eigen::MatrixXd::Identity(matrixToInvert.rows(), matrixToInvert.cols()));
  }
  else if(inversionMethod == "ldlt")
  {
    Eigen::LLT<Eigen::MatrixXd> lltOfA(matrixToInvert);
    return lltOfA.solve(Eigen::MatrixXd::Identity(matrixToInvert.rows(), matrixToInvert.cols()));
  }
  else if(inversionMethod == "bdcSvd")
  {
    return matrixToInvert.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(Eigen::MatrixXd::Identity(matrixToInvert.rows(), matrixToInvert.cols()));
  } 
  else{
    Info << "Unkwown inversion method, solving with : completeOrthogonalDecomposition" << endl;
    return matrixToInvert.completeOrthogonalDecomposition().pseudoInverse();
  }
}
 
 
}