Skip to content
Permalink
Browse files

Merge pull request #13405 from hugary1995/next

added AD compatible spectral decomposition for ranktwotensor
  • Loading branch information...
dschwen committed May 31, 2019
2 parents c3c0499 + 5b7b665 commit 5d1ad78212fb38b837fd75f7972a948bcb1e8fe1
@@ -369,6 +369,12 @@ class RankTwoTensorTempl : public TensorValue<T>
/// Print the rank two tensor
void print(std::ostream & stm = Moose::out) const;

/// Print the Real part of the DualReal rank two tensor
void printReal(std::ostream & stm = Moose::out) const;

/// Print the Real part of the DualReal rank two tensor along with its first nDual dual numbers
void printDualReal(unsigned int nDual, std::ostream & stm = Moose::out) const;

/// Add identity times a to _coords
void addIa(const T & a);

@@ -387,6 +393,43 @@ class RankTwoTensorTempl : public TensorValue<T>
*/
void symmetricEigenvalues(std::vector<T> & eigvals) const;

/**
* computes and returns the permutation matrix P
* @param old_elements is the original row/column numbering
* @param new_elements is the permuted row/column numbering
* Dual numbers are permuted as well
* P * A permutes rows and A * P^T permutes columns
*/
RankTwoTensorTempl<T>
permutationTensor(const std::array<unsigned int, LIBMESH_DIM> & old_elements,
const std::array<unsigned int, LIBMESH_DIM> & new_elements) const;

/**
* computes and returns the Givens rotation matrix R
* @param row1 is the row number of the first component to rotate
* @param row2 is the row number of the second component to rotate
* @param col is the column number of the components to rotate
* consider a RankTwoTensor A = [ a11 a12 a13
* a21 a22 a23
* a31 a32 a33]
* and we want to rotate a21 and a31. Then row1 = 1, row2 = 2, col = 0.
* It retunrs the Givens rotation matrix R of this tensor A such that R * A rotates the second
* component to zero, i.e. R * A = [ a11 a12 a13
* r a22 a23
* 0 a32 a33]
* A DualReal instantiation is available to rotate dual numbers as well.
*/
RankTwoTensorTempl<T>
givensRotation(unsigned int row1, unsigned int row2, unsigned int col) const;

/// computes the Hessenberg form of this matrix A and its unitary transformation U such that A = U * H * U^T
void hessenberg(RankTwoTensorTempl<T> & H, RankTwoTensorTempl<T> & U) const;

/// computes the QR factorization such that A = Q * R, where Q is the unitary matrix and R an upper triangular matrix
void QR(RankTwoTensorTempl<T> & Q,
RankTwoTensorTempl<T> & R,
unsigned int dim = RankTwoTensorTempl<T>::N) const;

/**
* computes eigenvalues and eigenvectors, assuming tens is symmetric, and places them
* in ascending order in eigvals. eigvecs is a matrix with the first column
@@ -829,6 +829,44 @@ RankTwoTensorTempl<T>::print(std::ostream & stm) const
}
}

template <>
void
RankTwoTensorTempl<Real>::printReal(std::ostream & stm) const
{
this->print(stm);
}

template <>
void
RankTwoTensorTempl<DualReal>::printReal(std::ostream & stm) const
{
const RankTwoTensorTempl<DualReal> & a = *this;
for (unsigned int i = 0; i < N; ++i)
{
for (unsigned int j = 0; j < N; ++j)
stm << std::setw(15) << a(i, j).value() << ' ';
stm << std::endl;
}
}

template <>
void
RankTwoTensorTempl<DualReal>::printDualReal(unsigned int nDual, std::ostream & stm) const
{
const RankTwoTensorTempl<DualReal> & a = *this;
for (unsigned int i = 0; i < N; ++i)
{
for (unsigned int j = 0; j < N; ++j)
{
stm << std::setw(15) << a(i, j).value() << " {";
for (unsigned int k = 0; k < nDual; ++k)
stm << std::setw(5) << a(i, j).derivatives()[k] << ' ';
stm << " }";
}
stm << std::endl;
}
}

template <typename T>
void
RankTwoTensorTempl<T>::addIa(const T & a)
@@ -889,6 +927,182 @@ RankTwoTensorTempl<T>::symmetricEigenvaluesEigenvectors(std::vector<T> & eigvals
eigvecs(j, i) = a[i * N + j];
}

template <typename T>
RankTwoTensorTempl<T>
RankTwoTensorTempl<T>::permutationTensor(
const std::array<unsigned int, LIBMESH_DIM> & old_elements,
const std::array<unsigned int, LIBMESH_DIM> & new_elements) const
{
RankTwoTensorTempl<T> P;

for (unsigned int i = 0; i < N; ++i)
P(old_elements[i], new_elements[i]) = 1.0;

return P;
}

template <typename T>
RankTwoTensorTempl<T>
RankTwoTensorTempl<T>::givensRotation(unsigned int row1, unsigned int row2, unsigned int col) const
{
T c, s;
T a = (*this)(row1, col);
T b = (*this)(row2, col);

if (MooseUtils::absoluteFuzzyEqual(b, 0.0) && MooseUtils::absoluteFuzzyEqual(a, 0.0))
{
c = a >= 0.0 ? 1.0 : -1.0;
s = 0.0;
}
else if (std::abs(a) > std::abs(b))
{
T t = b / a;
Real sgn = a >= 0.0 ? 1.0 : -1.0;
T u = sgn * std::sqrt(1.0 + t * t);
c = 1.0 / u;
s = c * t;
}
else
{
T t = a / b;
Real sgn = b >= 0.0 ? 1.0 : -1.0;
T u = sgn * std::sqrt(1.0 + t * t);
s = 1.0 / u;
c = s * t;
}

RankTwoTensorTempl<T> R(initIdentity);
R(row1, row1) = c;
R(row1, row2) = s;
R(row2, row1) = -s;
R(row2, row2) = c;

return R;
}

template <typename T>
void
RankTwoTensorTempl<T>::hessenberg(RankTwoTensorTempl<T> & H, RankTwoTensorTempl<T> & U) const
{
H = *this;
U.zero();
U.addIa(1.0);

if (N < 3)
return;

RankTwoTensorTempl<T> R = this->givensRotation(N - 2, N - 1, 0);
H = R * H * R.transpose();
U = U * R.transpose();
}

template <typename T>
void
RankTwoTensorTempl<T>::QR(RankTwoTensorTempl<T> & Q,
RankTwoTensorTempl<T> & R,
unsigned int dim) const
{
R = *this;
Q.zero();
Q.addIa(1.0);

for (unsigned int i = 0; i < dim - 1; i++)
for (unsigned int b = dim - 1; b > i; b--)
{
unsigned int a = b - 1;
RankTwoTensorTempl<T> CS = R.givensRotation(a, b, i);
R = CS * R;
Q = Q * CS.transpose();
}
}

template <>
void
RankTwoTensorTempl<DualReal>::QR(RankTwoTensorTempl<DualReal> & Q,
RankTwoTensorTempl<DualReal> & R,
unsigned int dim) const
{
R = *this;
Q.zero();
Q.addIa(1.0);

for (unsigned int i = 0; i < dim - 1; i++)
for (unsigned int b = dim - 1; b > i; b--)
{
unsigned int a = b - 1;

// special case when both entries to rotate are zero
// in which case the dual numbers cannot be rotated
// therefore we need to find another nonzero entry to permute
RankTwoTensorTempl<DualReal> P(initIdentity);
if (MooseUtils::absoluteFuzzyEqual(R(a, i).value(), 0.0) &&
MooseUtils::absoluteFuzzyEqual(R(b, i).value(), 0.0))
{
unsigned int c = 3 - a - b;
if (!MooseUtils::absoluteFuzzyEqual(R(c, i).value(), 0.0))
P = this->permutationTensor({a, b, c}, {c, b, a});
}

Q = Q * P.transpose();
R = P * R;
RankTwoTensorTempl<DualReal> CS = R.givensRotation(a, b, i);
R = P.transpose() * CS * R;
Q = Q * CS.transpose() * P;
}
}

template <>
void
RankTwoTensorTempl<DualReal>::symmetricEigenvaluesEigenvectors(
std::vector<DualReal> & eigvals, RankTwoTensorTempl<DualReal> & eigvecs) const
{
const Real eps = libMesh::TOLERANCE * libMesh::TOLERANCE;

eigvals.resize(N);
RankTwoTensorTempl<DualReal> D, Q, R;
this->hessenberg(D, eigvecs);

unsigned int iter = 0;
for (unsigned m = N - 1; m > 0; m--)
do
{
iter++;
DualReal shift = D(m, m);
D.addIa(-shift);
D.QR(Q, R, m + 1);
D = R * Q;
D.addIa(shift);
eigvecs = eigvecs * Q;
} while (std::abs(D(m, m - 1)) > eps);

for (unsigned int i = 0; i < N; i++)
eigvals[i] = D(i, i);

// Sort eigenvalues and corresponding vectors.
for (unsigned int i = 0; i < N - 1; i++)
{
unsigned int k = i;
DualReal p = eigvals[i];
for (unsigned int j = i + 1; j < N; j++)
if (eigvals[j] < p)
{
k = j;
p = eigvals[j];
}
if (k != i)
{
eigvals[k] = eigvals[i];
eigvals[i] = p;
for (unsigned int j = 0; j < N; j++)
{
p = eigvecs(j, i);
eigvecs(j, i) = eigvecs(j, k);
eigvecs(j, k) = p;
}
}
}
}

template <typename T>
void
RankTwoTensorTempl<T>::dsymmetricEigenvalues(std::vector<T> & eigvals,
@@ -982,8 +1196,8 @@ RankTwoTensorTempl<T>::syev(const char * calculation_type,

for (unsigned int i = 0; i < N; ++i)
for (unsigned int j = 0; j < N; ++j)
// a is destroyed by dsyev, and if calculation_type == "V" then eigenvectors are placed there
// Note the explicit symmeterisation
// a is destroyed by dsyev, and if calculation_type == "V" then eigenvectors are placed
// there Note the explicit symmeterisation
a[i * N + j] = 0.5 * (this->operator()(i, j) + this->operator()(j, i));

// compute the eigenvalues only (if calculation_type == "N"),

0 comments on commit 5d1ad78

Please sign in to comment.
You can’t perform that action at this time.