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 "SymmetricRankTwoTensorImplementation.h"
      11             : 
      12             : template <>
      13             : void
      14           2 : SymmetricRankTwoTensor::syev(const char * calculation_type,
      15             :                              std::vector<Real> & eigvals,
      16             :                              std::vector<Real> & a) const
      17             : {
      18           2 :   eigvals.resize(Ndim);
      19           2 :   a.resize(Ndim * Ndim);
      20             : 
      21             :   // prepare data for the LAPACKsyev_ routine (which comes from petscblaslapack.h)
      22           2 :   PetscBLASInt nd = Ndim;
      23           2 :   PetscBLASInt lwork = 66 * nd;
      24             :   PetscBLASInt info;
      25           2 :   std::vector<PetscScalar> work(lwork);
      26             : 
      27           2 :   auto A = RankTwoTensor(*this);
      28           8 :   for (auto i : make_range(Ndim))
      29          24 :     for (auto j : make_range(Ndim))
      30             :       // a is destroyed by dsyev, and if calculation_type == "V" then eigenvectors are placed
      31             :       // there
      32          18 :       a[i * Ndim + j] = A(i, j);
      33             : 
      34             :   // compute the eigenvalues only (if calculation_type == "N"),
      35             :   // or both the eigenvalues and eigenvectors (if calculation_type == "V")
      36             :   // assume upper triangle of a is stored (second "U")
      37           2 :   LAPACKsyev_(calculation_type, "U", &nd, &a[0], &nd, &eigvals[0], &work[0], &lwork, &info);
      38             : 
      39           2 :   if (info != 0)
      40           0 :     mooseError("In computing the eigenvalues and eigenvectors for the symmetric rank-2 tensor (",
      41           0 :                Moose::stringify(a),
      42             :                "), the PETSC LAPACK syev routine returned error code ",
      43             :                info);
      44           2 : }
      45             : 
      46             : template <>
      47             : void
      48           0 : SymmetricRankTwoTensor::symmetricEigenvalues(std::vector<Real> & eigvals) const
      49             : {
      50           0 :   std::vector<Real> a;
      51           0 :   syev("N", eigvals, a);
      52           0 : }
      53             : 
      54             : template <>
      55             : void
      56           2 : SymmetricRankTwoTensor::symmetricEigenvaluesEigenvectors(std::vector<Real> & eigvals,
      57             :                                                          RankTwoTensor & eigvecs) const
      58             : {
      59           2 :   std::vector<Real> a;
      60           2 :   syev("V", eigvals, a);
      61             : 
      62           8 :   for (auto i : make_range(Ndim))
      63          24 :     for (auto j : make_range(Ndim))
      64          18 :       eigvecs(j, i) = a[i * Ndim + j];
      65           2 : }