Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "RankTwoTensorImplementation.h"
      11             : 
      12             : template <>
      13             : void
      14           0 : RankTwoTensor::printReal(std::ostream & stm) const
      15             : {
      16           0 :   this->print(stm);
      17           0 : }
      18             : 
      19             : template <>
      20             : void
      21           0 : ADRankTwoTensor::printReal(std::ostream & stm) const
      22             : {
      23           0 :   const ADRankTwoTensor & a = *this;
      24           0 :   for (const auto i : make_range(N))
      25             :   {
      26           0 :     for (const auto j : make_range(N))
      27           0 :       stm << std::setw(15) << a(i, j).value() << ' ';
      28           0 :     stm << std::endl;
      29             :   }
      30           0 : }
      31             : 
      32             : template <>
      33             : void
      34           0 : ADRankTwoTensor::printADReal(unsigned int nDual, std::ostream & stm) const
      35             : {
      36           0 :   const ADRankTwoTensor & a = *this;
      37           0 :   for (const auto i : make_range(N))
      38             :   {
      39           0 :     for (const auto j : make_range(N))
      40             :     {
      41           0 :       stm << std::setw(15) << a(i, j).value() << " {";
      42           0 :       for (const auto k : make_range(nDual))
      43           0 :         stm << std::setw(5) << a(i, j).derivatives()[k] << ' ';
      44           0 :       stm << " }";
      45             :     }
      46           0 :     stm << std::endl;
      47             :   }
      48           0 : }
      49             : 
      50             : template <>
      51             : void
      52         630 : RankTwoTensor::syev(const char * calculation_type,
      53             :                     std::vector<Real> & eigvals,
      54             :                     std::vector<Real> & a) const
      55             : {
      56         630 :   eigvals.resize(N);
      57         630 :   a.resize(N * N);
      58             : 
      59             :   // prepare data for the LAPACKsyev_ routine (which comes from petscblaslapack.h)
      60         630 :   PetscBLASInt nd = N;
      61         630 :   PetscBLASInt lwork = 66 * nd;
      62             :   PetscBLASInt info;
      63         630 :   std::vector<PetscScalar> work(lwork);
      64             : 
      65        2520 :   for (const auto i : make_range(N))
      66        7560 :     for (const auto j : make_range(N))
      67             :       // a is destroyed by dsyev, and if calculation_type == "V" then eigenvectors are placed
      68             :       // there Note the explicit symmeterisation
      69        5670 :       a[i * N + j] = 0.5 * (this->operator()(i, j) + this->operator()(j, i));
      70             : 
      71             :   // compute the eigenvalues only (if calculation_type == "N"),
      72             :   // or both the eigenvalues and eigenvectors (if calculation_type == "V")
      73             :   // assume upper triangle of a is stored (second "U")
      74         630 :   LAPACKsyev_(calculation_type, "U", &nd, &a[0], &nd, &eigvals[0], &work[0], &lwork, &info);
      75             : 
      76         630 :   if (info != 0)
      77           0 :     mooseException(
      78             :         "In computing the eigenvalues and eigenvectors for the symmetric rank-2 tensor (",
      79             :         Moose::stringify(a),
      80             :         "), the PETSC LAPACK syev routine returned error code ",
      81             :         info);
      82         630 : }
      83             : 
      84             : template <>
      85             : void
      86          96 : RankTwoTensor::symmetricEigenvalues(std::vector<Real> & eigvals) const
      87             : {
      88          96 :   std::vector<Real> a;
      89          96 :   syev("N", eigvals, a);
      90          96 : }
      91             : 
      92             : template <>
      93             : void
      94          36 : RankTwoTensor::symmetricEigenvaluesEigenvectors(std::vector<Real> & eigvals,
      95             :                                                 RankTwoTensor & eigvecs) const
      96             : {
      97          36 :   std::vector<Real> a;
      98          36 :   syev("V", eigvals, a);
      99             : 
     100         144 :   for (const auto i : make_range(N))
     101         432 :     for (const auto j : make_range(N))
     102         324 :       eigvecs(j, i) = a[i * N + j];
     103          36 : }
     104             : 
     105             : template <>
     106             : void
     107          25 : ADRankTwoTensor::symmetricEigenvaluesEigenvectors(std::vector<ADReal> & eigvals,
     108             :                                                   ADRankTwoTensor & eigvecs) const
     109             : {
     110             :   typedef Eigen::Matrix<ADReal, N, N> RankTwoMatrix;
     111          25 :   RankTwoMatrix self;
     112         100 :   for (const auto i : make_range(N))
     113         225 :     for (unsigned int j = i; j < N; ++j)
     114             :     {
     115         150 :       auto & v = self(j, i);
     116         150 :       v = (*this)(i, j);
     117         150 :       if (i != j && MooseUtils::absoluteFuzzyEqual(v, 0.0))
     118          20 :         v.value() = 0.0;
     119             :     }
     120             : 
     121          25 :   Eigen::SelfAdjointEigenSolver<RankTwoMatrix> es;
     122          25 :   es.compute(self);
     123             : 
     124          25 :   const auto & lambda = es.eigenvalues();
     125          25 :   eigvals.resize(N);
     126         100 :   for (const auto i : make_range(N))
     127          75 :     eigvals[i] = lambda(i);
     128             : 
     129          25 :   const auto & v = es.eigenvectors();
     130         100 :   for (const auto i : make_range(N))
     131         300 :     for (const auto j : make_range(N))
     132         225 :       eigvecs(i, j) = v(i, j);
     133          25 : }
     134             : 
     135             : template <>
     136             : void
     137           0 : RankTwoTensor::getRUDecompositionRotation(RankTwoTensor & rot) const
     138             : {
     139           0 :   const RankTwoTensor & a = *this;
     140           0 :   RankTwoTensor c, diag, evec;
     141             :   PetscScalar cmat[N][N], work[10];
     142             :   PetscReal w[N];
     143             : 
     144             :   // prepare data for the LAPACKsyev_ routine (which comes from petscblaslapack.h)
     145           0 :   PetscBLASInt nd = N, lwork = 10, info;
     146             : 
     147           0 :   c = a.transpose() * a;
     148             : 
     149           0 :   for (const auto i : make_range(N))
     150           0 :     for (const auto j : make_range(N))
     151           0 :       cmat[i][j] = c(i, j);
     152             : 
     153           0 :   LAPACKsyev_("V", "U", &nd, &cmat[0][0], &nd, w, work, &lwork, &info);
     154             : 
     155           0 :   if (info != 0)
     156           0 :     mooseException(
     157             :         "In computing the eigenvalues and eigenvectors of a symmetric rank-2 tensor, the "
     158             :         "PETSC LAPACK syev routine returned error code ",
     159             :         info);
     160             : 
     161           0 :   diag.zero();
     162             : 
     163           0 :   for (const auto i : make_range(N))
     164           0 :     diag(i, i) = std::sqrt(w[i]);
     165             : 
     166           0 :   for (const auto i : make_range(N))
     167           0 :     for (const auto j : make_range(N))
     168           0 :       evec(i, j) = cmat[i][j];
     169             : 
     170           0 :   rot = a * (evec.transpose() * diag * evec).inverse();
     171           0 : }