From: Jean-Paul Pelteret Date: Fri, 28 Jul 2017 16:40:34 +0000 (+0200) Subject: Split header and explicitly instantiate template functions. X-Git-Tag: v9.0.0-rc1~1366^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=23383eddd5953628b50232195f7e14e4f6df034f;p=dealii.git Split header and explicitly instantiate template functions. --- diff --git a/include/deal.II/base/symmetric_tensor.h b/include/deal.II/base/symmetric_tensor.h index f4bc00c5a6..2ded82e272 100644 --- a/include/deal.II/base/symmetric_tensor.h +++ b/include/deal.II/base/symmetric_tensor.h @@ -2326,7 +2326,6 @@ Number second_invariant (const SymmetricTensor<2,3,Number> &t) /** * Return the eigenvalues of a symmetric 1x1 tensor of rank 2. - * * The (single) entry of the tensor is, of course, equal to the (single) * eigenvalue. * @@ -2335,10 +2334,7 @@ Number second_invariant (const SymmetricTensor<2,3,Number> &t) */ template std::array -eigenvalues (const SymmetricTensor<2,1,Number> &T) -{ - return { {T[0][0]} }; -} +eigenvalues (const SymmetricTensor<2,1,Number> &T); @@ -2365,40 +2361,7 @@ eigenvalues (const SymmetricTensor<2,1,Number> &T) */ template std::array -eigenvalues (const SymmetricTensor<2,2,Number> &T) -{ - const Number upp_tri_sq = T[0][1]*T[0][1]; - if (upp_tri_sq == Number(0.0)) - { - // The tensor is diagonal - std::array eig_vals = - { - {T[0][0], T[1][1]} - }; - - // Sort from largest to smallest. - std::sort(eig_vals.begin(), eig_vals.end(), std::greater()); - return eig_vals; - } - else - { - const Number tr_T = trace(T); - const Number det_T = determinant(T); - const Number descrim = tr_T*tr_T - 4.0*det_T; - Assert(descrim > Number(0.0), ExcMessage("The roots of the characteristic polynomial are complex valued.")); - const Number sqrt_desc = std::sqrt(descrim); - - std::array eig_vals = - { - { - 0.5*(tr_T + sqrt_desc), - 0.5*(tr_T - sqrt_desc) - } - }; - Assert(eig_vals[0] >= eig_vals[1], ExcMessage("The eigenvalue ordering is incorrect.")); - return eig_vals; - } -} +eigenvalues (const SymmetricTensor<2,2,Number> &T); @@ -2423,511 +2386,154 @@ eigenvalues (const SymmetricTensor<2,2,Number> &T) */ template std::array -eigenvalues (const SymmetricTensor<2,3,Number> &T) -{ - const Number upp_tri_sq = T[0][1]*T[0][1] + T[0][2]*T[0][2] + T[1][2]*T[1][2]; - if (upp_tri_sq == Number(0.0)) - { - // The tensor is diagonal - std::array eig_vals - = { {T[0][0], T[1][1], T[2][2]} }; - - // Sort from largest to smallest. - std::sort(eig_vals.begin(), eig_vals.end(), std::greater()); - return eig_vals; - } - else - { - // Perform an affine change to T, and solve a different - // characteristic equation that has a trigonometric solution. - // Decompose T = p*B + q*I , and set q = tr(T)/3 - // and p = (tr((T - q.I)^{2})/6)^{1/2} . Then solve the equation - // 0 = det(\lambda*I - B) = \lambda^{3} - 3*\lambda - det(B) - // which has the solution - // \lambda = 2*cos(1/3 * acos(det(B)/2) +2/3*pi*k ) ; k = 0,1,2 - // when substituting \lambda = 2.cos(theta) and using trig identities. - const Number tr_T = trace(T); - const Number q = tr_T/3.0; - const Number tmp1 = ( T[0][0] - q)*(T[0][0] - q) - + (T[1][1] - q)*(T[1][1] - q) - + (T[2][2] - q)*(T[2][2] - q) - + 2.0 * upp_tri_sq; - const Number p = std::sqrt(tmp1/6.0); - const SymmetricTensor<2,3,Number> B = (1.0/p)*(T - q*unit_symmetric_tensor<3,Number>()); - const Number tmp_2 = determinant(B)/2.0; - - // The value of tmp_2 should be within [-1,1], however - // floating point errors might place it slightly outside - // this range. It is therefore necessary to correct for it - const Number phi = - (tmp_2 <= -1.0 ? M_PI/3.0 : - (tmp_2 >= 1.0 ? 0.0 : - std::acos(tmp_2)/3.0)); - - // Due to the trigonometric solution, the computed eigenvalues - // should be predictably in the order eig1 >= eig2 >= eig3... - std::array eig_vals - = { { - q + 2.0*p *std::cos(phi), - 0.0, - q + 2.0*p *std::cos(phi + (2.0/3.0*M_PI)) - } - }; - // Use the identity tr(T) = eig1 + eig2 + eig3 - eig_vals[1] = tr_T - eig_vals[0] - eig_vals[2]; - - // ... however, when equal roots exist then floating point - // errors may make this no longer be the case. - // Sort from largest to smallest. - std::sort(eig_vals.begin(), eig_vals.end(), std::greater()); - - return eig_vals; - } -} +eigenvalues (const SymmetricTensor<2,3,Number> &T); namespace internal { - namespace + /** + * A namespace for functions and classes that are internal to how the + * SymmetricTensor class (and its associate functions) works. + */ + namespace SymmetricTensor { /** - * Tridiagonalize a rank-2 symmetric using the Householder method. + * Tridiagonalize a rank-2 symmetric tensor using the Householder method. * The specialized algorithm implemented here is given in - * Kopp, J. - * Efficient numerical diagonalization of hermitian 3x3 matrices - * International Journal of Modern Physics C, 2008, 19, 523-548 - * doi: 10.1142/S0129183108012303 - * arXiv.org preprint: physics/0610206 + * @code{.bib} + * @Article{Kopp2008, + * title = {Efficient numerical diagonalization of hermitian 3x3 matrices}, + * author = {Kopp, J.}, + * journal = {International Journal of Modern Physics C}, + * year = {2008}, + * volume = {19}, + * number = {3}, + * pages = {523--548}, + * doi = {10.1142/S0129183108012303}, + * eprinttype = {arXiv}, + * eprint = {physics/0610206v3}, + * eprintclass = {physics.comp-ph}, + * url = {https://www.mpi-hd.mpg.de/personalhomes/globes/3x3/index.html} + * } + * @endcode * and is based off of the generic algorithm presented in section 11.3.2 of - * Press, W. H. - * Numerical recipes 3rd edition: The art of scientific computing - * Cambridge university press, 2007 + * @code{.bib} + * @Book{Press2007, + * title = {Numerical recipes 3rd edition: The art of scientific computing}, + * author = {Press, W. H.}, + * journal = {Cambridge university press}, + * year = {2007} + * } + * @endcode * * @param[in] A This tensor to be tridiagonalized * @param[out] Q The orthogonal matrix effecting the transformation * @param[out] d The diagonal elements of the tridiagonal matrix * @param[out] e The off-diagonal elements of the tridiagonal matrix + * + * @author Joachim Kopp, Jean-Paul Pelteret, 2017 */ template void tridiagonalize (const dealii::SymmetricTensor<2,dim,Number> &A, dealii::Tensor<2,dim,Number> &Q, std::array &d, - std::array &e) - { - // Create some intermediate storage - Number h,g,omega_inv,K,f; - - // Initialize the transformation matrix as the - // identity tensor - Q = dealii::unit_symmetric_tensor(); - - // Make the first row and column to be of the - // desired form - h = 0.0; - for (int i=1; i < dim; i++) - h += A[0][i]*A[0][i]; - - g = 0.0; - if (A[0][1] > 0.0) - g = -std::sqrt(h); - else - g = std::sqrt(h); - e[0] = g; + std::array &e); - std::array u; - for (int i=1; i < dim; i++) - { - u[i] = A[0][i]; - if (i == 1) - u[i] -= g; - } - - std::array q; - const Number omega = h - g * A[0][1]; - if (omega > 0.0) - { - omega_inv = 1.0 / omega; - K = 0.0; - for (int i=1; i < dim; i++) - { - f = 0.0; - for (int j=1; j < dim; j++) - f += A[i][j] * u[j]; - q[i] = omega_inv * f; - K += u[i] * f; - } - K *= 0.5*omega_inv*omega_inv; - - for (int i=1; i < dim; i++) - q[i] = q[i] - K * u[i]; - - d[0] = A[0][0]; - for (int i=1; i < dim; i++) - d[i] = A[i][i] - 2.0*q[i]*u[i]; - - // Store inverse Householder transformation - // in Q - for (int j=1; j < dim; j++) - { - f = omega_inv * u[j]; - for (int i=1; i < dim; i++) - Q[i][j] = Q[i][j] - f*u[i]; - } - - // For dim = 3: Calculate updated A[1][2] and - // store it in e[1] - for (int i=1; i < dim-1; i++) - e[i] = A[i][i+1] - q[i]*u[i+1] - u[i]*q[i+1]; - } - else - { - for (int i=0; i < dim; i++) - d[i] = A[i][i]; - - // For dim = 3: - for (int i=1; i < dim-1; i++) - e[i] = A[i][i+1]; - } - } /** * Compute the eigenvalues and eigenvectors of a real-valued rank-2 * symmetric tensor using the QL algorithm with implicit shifts. * The specialized algorithm implemented here is given in - * Kopp, J. - * Efficient numerical diagonalization of hermitian 3x3 matrices - * International Journal of Modern Physics C, 2008, 19, 523-548 - * doi: 10.1142/S0129183108012303 - * arXiv.org preprint: physics/0610206 + * @code{.bib} + * @Article{Kopp2008, + * title = {Efficient numerical diagonalization of hermitian 3x3 matrices}, + * author = {Kopp, J.}, + * journal = {International Journal of Modern Physics C}, + * year = {2008}, + * volume = {19}, + * number = {3}, + * pages = {523--548}, + * doi = {10.1142/S0129183108012303}, + * eprinttype = {arXiv}, + * eprint = {physics/0610206v3}, + * eprintclass = {physics.comp-ph}, + * url = {https://www.mpi-hd.mpg.de/personalhomes/globes/3x3/index.html} + * } + * @endcode * and is based off of the generic algorithm presented in section 11.4.3 of - * Press, W. H. - * Numerical recipes 3rd edition: The art of scientific computing - * Cambridge university press, 2007. + * @code{.bib} + * @Book{Press2007, + * title = {Numerical recipes 3rd edition: The art of scientific computing}, + * author = {Press, W. H.}, + * journal = {Cambridge university press}, + * year = {2007} + * } + * @endcode * * @param[in] A The tensor of which the eigenvectors and eigenvalues are * to be computed. * - * @return An array containing the eigenvectors and the associated eigenvalues + * @return An array containing the eigenvectors and the associated eigenvalues. + * The array is not sorted in any particular order. + * + * @author Joachim Kopp, Jean-Paul Pelteret, 2017 */ template std::array >,dim> - ql_implicit_shifts (const dealii::SymmetricTensor<2,dim,Number> &A) - { - static_assert(numbers::NumberTraits::is_complex == false, - "This implementation of the QL implicit shift algorithm does " - "not support complex numbers"); - - // Transform A to real tridiagonal form by the Householder method: - // The orthogonal matrix effecting the transformation - // this will ultimately store the eigenvectors - dealii::Tensor<2,dim,Number> Q; - // The diagonal elements of the tridiagonal matrix; - // this will ultimately store the eigenvalues - std::array w; - // The off-diagonal elements of the tridiagonal - std::array ee; - tridiagonalize(A, Q, w, ee); - - // Number of iterations - const unsigned int max_n_it = 30; - - // Transfer the off-diagonal entries to an auxiliary array - // The third element is used only as temporary workspace - std::array e; - for (unsigned int i=0; i >,dim> (); - } - - // Calculate the shift.. - g = (w[l+1] - w[l]) / (e[l] + e[l]); - r = std::sqrt(g*g + 1.0); - // .. and then compute g = d_m - k_s for the - // plane rotation (Press2007a eq 11.4.22) - if (g > 0.0) - g = w[m] - w[l] + e[l]/(g + r); - else - g = w[m] - w[l] + e[l]/(g - r); - - // Perform plane rotation, as is done in the - // standard QL algorithm, followed by Givens - // rotations to recover the tridiagonal form - s = c = 1.0; - p = 0.0; - for (int i=m-1; i >= l; i--) - { - f = s * e[i]; - b = c * e[i]; - - // Branch to recover from underflow - if (std::abs(f) > std::abs(g)) - { - c = g / f; - r = std::sqrt(c*c + 1.0); - e[i+1] = f * r; - c *= (s = 1.0/r); - } - else - { - s = f / g; - r = std::sqrt(s*s + 1.0); - e[i+1] = g * r; - s *= (c = 1.0/r); - } - - g = w[i+1] - p; - r = (w[i] - g)*s + 2.0*c*b; - p = s * r; - w[i+1] = g + p; - g = c*r - b; - - // Form the eigenvectors - for (int k=0; k < dim; k++) - { - t = Q[k][i+1]; - Q[k][i+1] = s*Q[k][i] + c*t; - Q[k][i] = c*Q[k][i] - s*t; - } - } - w[l] -= p; - e[l] = g; - e[m] = 0.0; - } - } + ql_implicit_shifts (const dealii::SymmetricTensor<2,dim,Number> &A); - // Structure the data to be outputted - std::array >,dim> eig_vals_vecs; - for (unsigned int e=0; e std::array >,dim> - jacobi (dealii::SymmetricTensor<2,dim,Number> A) - { - static_assert(numbers::NumberTraits::is_complex == false, - "This implementation of the Jacobi algorithm does " - "not support complex numbers"); - - // Sums of diagonal resp. off-diagonal elements - Number sd, so; - // sin(phi), cos(phi), tan(phi) and temporary storage - Number s, c, t; - // More temporary storage - Number g, h, z, theta; - // Threshold value - Number thresh; - - // Initialize the transformation matrix as the - // identity tensor - dealii::Tensor<2,dim,Number> Q (dealii::unit_symmetric_tensor()); - - // The diagonal elements of the tridiagonal matrix; - // this will ultimately store the eigenvalues - std::array w; - for (int i=0; i < dim; i++) - w[i] = A[i][i]; - - // Calculate (tr(A))^{2} - sd = trace(A); - sd *= sd; - - // Number of iterations - const unsigned int max_n_it = 50; - for (unsigned int it=0; it <= max_n_it; it++) - { - // Test for convergence - so = 0.0; - for (int p=0; p < dim; p++) - for (int q=p+1; q < dim; q++) - so += std::abs(A[p][q]); - if (so == 0.0) - break; - - // Throw if no convergence is achieved within a - // stipulated number of iterations - if (it == max_n_it) - { - AssertThrow(false, ExcMessage("No convergence in iterative Jacobi eigenvector algorithm.")) - return std::array >,dim> (); - } - - // Compute threshold value which dictates whether or - // not a Jacobi rotation is performed - const unsigned int n_it_skip = 4; - if (it < n_it_skip) - thresh = 0.2 * so / (dim*dim); - else - thresh = 0.0; - - // Perform sweep - for (int p=0; p < dim; p++) - for (int q=p+1; q < dim; q++) - { - g = 100.0 * std::abs(A[p][q]); - - // After a given number of iterations the - // rotation is skipped if the off-diagonal - // element is small - if (it > n_it_skip && - std::abs(w[p]) + g == std::abs(w[p]) && - std::abs(w[q]) + g == std::abs(w[q])) - { - A[p][q] = 0.0; - } - else if (std::abs(A[p][q]) > thresh) - { - // Calculate Jacobi transformation - h = w[q] - w[p]; - - // Compute surrogate for angle theta resulting from - // angle transformation and subsequent smallest solution - // of quadratic equation - if (std::abs(h) + g == std::abs(h)) - { - // Prevent overflow for large theta^2. This computation - // is the algebraic equivalent of t = 1/(2*theta). - t = A[p][q] / h; - } - else - { - theta = 0.5 * h / A[p][q]; - if (theta < 0.0) - t = -1.0 / (std::sqrt(1.0 + theta*theta) - theta); - else - t = 1.0 / (std::sqrt(1.0 + theta*theta) + theta); - } - - // Compute trigonometric functions for rotation - // in such a way as to prevent overflow for - // large theta. - c = 1.0/std::sqrt(1.0 + t*t); - s = t * c; - z = t * A[p][q]; - - // Apply Jacobi transformation... - A[p][q] = 0.0; - w[p] -= z; - w[q] += z; - // ... by executing the various rotations in sequence - for (int r=0; r < p; r++) - { - t = A[r][p]; - A[r][p] = c*t - s*A[r][q]; - A[r][q] = s*t + c*A[r][q]; - } - for (int r=p+1; r < q; r++) - { - t = A[p][r]; - A[p][r] = c*t - s*A[r][q]; - A[r][q] = s*t + c*A[r][q]; - } - for (int r=q+1; r < dim; r++) - { - t = A[p][r]; - A[p][r] = c*t - s*A[q][r]; - A[q][r] = s*t + c*A[q][r]; - } - - // Update the eigenvectors - for (int r=0; r < dim; r++) - { - t = Q[r][p]; - Q[r][p] = c*t - s*Q[r][q]; - Q[r][q] = s*t + c*Q[r][q]; - } - } - } - } - - // Structure the data to be outputted - std::array >,dim> eig_vals_vecs; - for (unsigned int e=0; e A); /** * Compute the eigenvalues and eigenvectors of a real-valued rank-2 - * symmetric tensor using the characteristic equation to compute eigenvalues + * symmetric 2x2 tensor using the characteristic equation to compute eigenvalues * and an analytical approach based on the cross-product for the eigenvectors. * If the computations are deemed too inaccurate then the method falls * back to ql_implicit_shifts. @@ -2935,228 +2541,88 @@ namespace internal * @param[in] A The tensor of which the eigenvectors and eigenvalues are * to be computed. * - * @return An array containing the eigenvectors and the associated eigenvalues + * @return An array containing the eigenvectors and the associated eigenvalues. + * The array is not sorted in any particular order. + * + * @author Joachim Kopp, Jean-Paul Pelteret, 2017 */ template std::array >,2> - hybrid (const dealii::SymmetricTensor<2,2,Number> &A) - { - static_assert(numbers::NumberTraits::is_complex == false, - "This implementation of the 2d Hybrid algorithm does " - "not support complex numbers"); - - const unsigned int dim = 2; - - // Calculate eigenvalues - const std::array w = eigenvalues(A); - - std::array >,dim> eig_vals_vecs; - - Number t, u; // Intermediate storage - t = std::abs(w[0]); - for (unsigned int i=1; i t) - t = u; - } - - if (t < 1.0) - u = t; - else - u = t*t; - - // Estimated maximum roundoff error - const double error = 256.0 * std::numeric_limits::epsilon() * u*u; - - // Store eigenvalues - eig_vals_vecs[0].first = w[0]; - eig_vals_vecs[1].first = w[1]; - - // Compute eigenvectors - // http://www.math.harvard.edu/archive/21b_fall_04/exhibits/2dmatrices/ - // https://math.stackexchange.com/a/1548616 - if (A[1][0] != 0.0) - { - // First eigenvector - eig_vals_vecs[0].second[0] = w[0] - A[1][1]; - eig_vals_vecs[0].second[1] = A[1][0]; - - // Second eigenvector - eig_vals_vecs[1].second[0] = w[1] - A[1][1]; - eig_vals_vecs[1].second[1] = A[1][0]; - } - else - { - // First eigenvector - eig_vals_vecs[0].second[0] = w[0]; - eig_vals_vecs[0].second[1] = 0.0; - - // Second eigenvector - eig_vals_vecs[1].second[0] = 0.0; - eig_vals_vecs[1].second[1] = w[1]; - } - // Normalize - eig_vals_vecs[0].second /= eig_vals_vecs[0].second.norm(); - eig_vals_vecs[1].second /= eig_vals_vecs[1].second.norm(); - - // If vectors are nearly linearly dependent, or if there might have - // been large cancelations in the calculation of A[i][i] - w[0], fall - // back to QL algorithm - if (eig_vals_vecs[0].second * eig_vals_vecs[1].second > error) - { - return ql_implicit_shifts(A); - } - - return eig_vals_vecs; - } + hybrid (const dealii::SymmetricTensor<2,2,Number> &A); /** * Compute the eigenvalues and eigenvectors of a real-valued rank-2 - * symmetric tensor using the characteristic equation to compute eigenvalues + * symmetric 3x3 tensor using the characteristic equation to compute eigenvalues * and an analytical approach based on the cross-product for the eigenvectors. * If the computations are deemed too inaccurate then the method falls * back to ql_implicit_shifts. * The specialized algorithm implemented here is given in - * Kopp, J. - * Efficient numerical diagonalization of hermitian 3x3 matrices - * International Journal of Modern Physics C, 2008, 19, 523-548 - * doi: 10.1142/S0129183108012303 - * arXiv.org preprint: physics/0610206 + * @code{.bib} + * @Article{Kopp2008, + * title = {Efficient numerical diagonalization of hermitian 3x3 matrices}, + * author = {Kopp, J.}, + * journal = {International Journal of Modern Physics C}, + * year = {2008}, + * volume = {19}, + * number = {3}, + * pages = {523--548}, + * doi = {10.1142/S0129183108012303}, + * eprinttype = {arXiv}, + * eprint = {physics/0610206v3}, + * eprintclass = {physics.comp-ph}, + * url = {https://www.mpi-hd.mpg.de/personalhomes/globes/3x3/index.html} + * } + * @endcode * * @param[in] A The tensor of which the eigenvectors and eigenvalues are * to be computed. * - * @return An array containing the eigenvectors and the associated eigenvalues + * @return An array containing the eigenvectors and the associated eigenvalues. + * The array is not sorted in any particular order. + * + * @author Joachim Kopp, Jean-Paul Pelteret, 2017 */ template std::array >,3> - hybrid (const dealii::SymmetricTensor<2,3,Number> &A) - { - static_assert(numbers::NumberTraits::is_complex == false, - "This implementation of the 3d Hybrid algorithm does " - "not support complex numbers"); + hybrid (const dealii::SymmetricTensor<2,3,Number> &A); - const unsigned int dim = 3; - Number norm; // Squared norm or inverse norm of current eigenvector - Number t, u; // Intermediate storage - - // Calculate eigenvalues - const std::array w = eigenvalues(A); - - t = std::abs(w[0]); - for (unsigned int i=1; i t) - t = u; - } - - if (t < 1.0) - u = t; - else - u = t*t; - - // Estimated maximum roundoff error - const double error = 256.0 * std::numeric_limits::epsilon() * u*u; - - // Initialize the transformation matrix as the - // identity tensor - dealii::Tensor<2,dim,Number> Q; - Q[0][1] = A[0][1]*A[1][2] - A[0][2]*A[1][1]; - Q[1][1] = A[0][2]*A[0][1] - A[1][2]*A[0][0]; - Q[2][1] = A[0][1]*A[0][1]; - - // Calculate first eigenvector by the formula - // v[0] = (A - w[0]).e1 x (A - w[0]).e2 - Q[0][0] = Q[0][1] + A[0][2]*w[0]; - Q[1][0] = Q[1][1] + A[1][2]*w[0]; - Q[2][0] = (A[0][0] - w[0]) * (A[1][1] - w[0]) - Q[2][1]; - norm = Q[0][0]*Q[0][0] + Q[1][0]*Q[1][0] + Q[2][0]*Q[2][0]; - - // If vectors are nearly linearly dependent, or if there might have - // been large cancellations in the calculation of A[i][i] - w[0], fall - // back to QL algorithm - // Note that this simultaneously ensures that multiple eigenvalues do - // not cause problems: If w[0] = w[1], then A - w[0] * I has rank 1, - // i.e. all columns of A - w[0] * I are linearly dependent. - if (norm <= error) - { - return ql_implicit_shifts(A); - } - else // This is the standard branch - { - norm = std::sqrt(1.0 / norm); - for (unsigned j=0; j < dim; j++) - Q[j][0] = Q[j][0] * norm; - } + namespace + { - // Calculate second eigenvector by the formula - // v[1] = (A - w[1]).e1 x (A - w[1]).e2 - Q[0][1] = Q[0][1] + A[0][2]*w[1]; - Q[1][1] = Q[1][1] + A[1][2]*w[1]; - Q[2][1] = (A[0][0] - w[1]) * (A[1][1] - w[1]) - Q[2][1]; - norm = Q[0][1]*Q[0][1] + Q[1][1]*Q[1][1] + Q[2][1]*Q[2][1]; - if (norm <= error) - { - return ql_implicit_shifts(A); - } - else + /** + * A struct that is used to sort arrays of pairs of eign=envalues and + * eigenvectors. Sorting is performed in in descending order of eigenvalue. + */ + template + struct SortEigenValuesVectors + { + typedef std::pair > EigValsVecs; + bool operator() (const EigValsVecs &lhs, + const EigValsVecs &rhs) { - norm = std::sqrt(1.0 / norm); - for (unsigned int j=0; j < dim; j++) - Q[j][1] = Q[j][1] * norm; + return lhs.first > rhs.first; } + }; - // Calculate third eigenvector according to - // v[2] = v[0] x v[1] - Q[0][2] = Q[1][0]*Q[2][1] - Q[2][0]*Q[1][1]; - Q[1][2] = Q[2][0]*Q[0][1] - Q[0][0]*Q[2][1]; - Q[2][2] = Q[0][0]*Q[1][1] - Q[1][0]*Q[0][1]; - - // Structure the data to be outputted - std::array >,dim> eig_vals_vecs; - for (unsigned int e=0; e - struct SortEigenValuesVectors - { - typedef std::pair > EigValsVecs; - bool operator() (const EigValsVecs &lhs, - const EigValsVecs &rhs) - { - return lhs.first > rhs.first; - } - }; - - } } // namespace internal +// The line below is to ensure that doxygen puts the full description +// of this global enumeration into the documentation +// See https://stackoverflow.com/a/1717984 +/** @file */ /** * An enumeration for the algorithm to be employed when performing * the computation of normalized eigenvectors and their corresponding - * eigenvalues. + * eigenvalues by the eigenvalues() and eigenvectors() methods operating on + * SymmetricTensor objects. * * The specialized algorithms utilized in computing the eigenvectors are * presented in @@ -3177,7 +2643,7 @@ namespace internal * } * @endcode */ -enum EigenvectorMethod +enum struct SymmetricTensorEigenvectorMethod { /** * A hybrid approach that preferentially uses the characteristic equation to @@ -3212,15 +2678,15 @@ enum EigenvectorMethod /** - * Return the eigenvalues and eigenvectors of a symmetric tensor of rank 2. + * Return the eigenvalues and eigenvectors of a symmetric 1x1 tensor of rank 2. * * @relates SymmetricTensor * @author Jean-Paul Pelteret, 2017 */ template std::array >,1> -eigenvectors (const SymmetricTensor<2,1,Number> &T, - const enum EigenvectorMethod /*method*/) +eigenvectors (const SymmetricTensor<2,1,Number> &T, + const enum SymmetricTensorEigenvectorMethod /*method*/ = SymmetricTensorEigenvectorMethod::ql_implicit_shifts) { return { {std::make_pair(T[0][0], Tensor<1,1,Number>({1.0}))} }; } @@ -3229,7 +2695,8 @@ eigenvectors (const SymmetricTensor<2,1,Number> &T, /** * Return the eigenvalues and eigenvectors of a real-valued rank-2 symmetric - * tensor $T$. + * tensor $T$. The array of matched eigenvalue and eigenvector pairs is sorted + * in descending order (determined by the eigenvalues). * * The specialized algorithms utilized in computing the eigenvectors are * presented in @@ -3253,73 +2720,31 @@ eigenvectors (const SymmetricTensor<2,1,Number> &T, * @relates SymmetricTensor * @author Joachim Kopp, Jean-Paul Pelteret, 2017 */ -template -std::array >,2> -eigenvectors (const SymmetricTensor<2,2,Number> &T, - const enum EigenvectorMethod method = ql_implicit_shifts) -{ - std::array >,2> eig_vals_vecs; - - if (method == hybrid) - eig_vals_vecs = internal::hybrid(T); - else if (method == ql_implicit_shifts) - eig_vals_vecs = internal::ql_implicit_shifts(T); - else if (method == jacobi) - eig_vals_vecs = internal::jacobi(T); - else - AssertThrow(false, ExcNotImplemented()); - - std::sort(eig_vals_vecs.begin(), eig_vals_vecs.end(), - internal::SortEigenValuesVectors<2,Number>()); - return eig_vals_vecs; -} - - - -/** - * Return the eigenvalues and eigenvectors of a real-valued rank-2 symmetric - * tensor $T$. - * - * The specialized algorithms utilized in computing the eigenvectors are - * presented in - * @code{.bib} - * @Article{Kopp2008, - * title = {Efficient numerical diagonalization of hermitian 3x3 matrices}, - * author = {Kopp, J.}, - * journal = {International Journal of Modern Physics C}, - * year = {2008}, - * volume = {19}, - * number = {3}, - * pages = {523--548}, - * doi = {10.1142/S0129183108012303}, - * eprinttype = {arXiv}, - * eprint = {physics/0610206v3}, - * eprintclass = {physics.comp-ph}, - * url = {https://www.mpi-hd.mpg.de/personalhomes/globes/3x3/index.html} - * } - * @endcode - * - * @relates SymmetricTensor - * @author Joachim Kopp, Jean-Paul Pelteret, 2017 - */ -template -std::array >,3> -eigenvectors (const SymmetricTensor<2,3,Number> &T, - const enum EigenvectorMethod method = ql_implicit_shifts) +template +std::array >,dim> +eigenvectors (const SymmetricTensor<2,dim,Number> &T, + const enum SymmetricTensorEigenvectorMethod method = SymmetricTensorEigenvectorMethod::ql_implicit_shifts) { - std::array >,3> eig_vals_vecs; + std::array >,dim> eig_vals_vecs; - if (method == hybrid) - eig_vals_vecs = internal::hybrid(T); - else if (method == ql_implicit_shifts) - eig_vals_vecs = internal::ql_implicit_shifts(T); - else if (method == jacobi) - eig_vals_vecs = internal::jacobi(T); - else - AssertThrow(false, ExcNotImplemented()); + switch (method) + { + case SymmetricTensorEigenvectorMethod::hybrid: + eig_vals_vecs = internal::SymmetricTensor::hybrid(T); + break; + case SymmetricTensorEigenvectorMethod::ql_implicit_shifts: + eig_vals_vecs = internal::SymmetricTensor::ql_implicit_shifts(T); + break; + case SymmetricTensorEigenvectorMethod::jacobi: + eig_vals_vecs = internal::SymmetricTensor::jacobi(T); + break; + default: + AssertThrow(false, ExcNotImplemented()); + } + // Sort in descending order before output. std::sort(eig_vals_vecs.begin(), eig_vals_vecs.end(), - internal::SortEigenValuesVectors<3,Number>()); + internal::SymmetricTensor::SortEigenValuesVectors()); return eig_vals_vecs; } diff --git a/include/deal.II/base/symmetric_tensor.templates.h b/include/deal.II/base/symmetric_tensor.templates.h new file mode 100644 index 0000000000..6c9e7fd3f9 --- /dev/null +++ b/include/deal.II/base/symmetric_tensor.templates.h @@ -0,0 +1,714 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, 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 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#ifndef dealii__symmetric_tensor_templates_h +#define dealii__symmetric_tensor_templates_h + + +#include + +#include + +DEAL_II_NAMESPACE_OPEN + + + +template +std::array +eigenvalues (const SymmetricTensor<2,1,Number> &T) +{ + return { {T[0][0]} }; +} + + + +template +std::array +eigenvalues (const SymmetricTensor<2,2,Number> &T) +{ + const Number upp_tri_sq = T[0][1]*T[0][1]; + if (upp_tri_sq == Number(0.0)) + { + // The tensor is diagonal + std::array eig_vals = + { + {T[0][0], T[1][1]} + }; + + // Sort from largest to smallest. + std::sort(eig_vals.begin(), eig_vals.end(), std::greater()); + return eig_vals; + } + else + { + const Number tr_T = trace(T); + const Number det_T = determinant(T); + const Number descrim = tr_T*tr_T - 4.0*det_T; + Assert(descrim > Number(0.0), ExcMessage("The roots of the characteristic polynomial are complex valued.")); + const Number sqrt_desc = std::sqrt(descrim); + + const std::array eig_vals = + { + { + static_cast(0.5*(tr_T + sqrt_desc)), + static_cast(0.5*(tr_T - sqrt_desc)) + } + }; + Assert(eig_vals[0] >= eig_vals[1], ExcMessage("The eigenvalue ordering is incorrect.")); + return eig_vals; + } +} + + + +template +std::array +eigenvalues (const SymmetricTensor<2,3,Number> &T) +{ + const Number upp_tri_sq = T[0][1]*T[0][1] + T[0][2]*T[0][2] + T[1][2]*T[1][2]; + if (upp_tri_sq == Number(0.0)) + { + // The tensor is diagonal + std::array eig_vals + = { {T[0][0], T[1][1], T[2][2]} }; + + // Sort from largest to smallest. + std::sort(eig_vals.begin(), eig_vals.end(), std::greater()); + return eig_vals; + } + else + { + // Perform an affine change to T, and solve a different + // characteristic equation that has a trigonometric solution. + // Decompose T = p*B + q*I , and set q = tr(T)/3 + // and p = (tr((T - q.I)^{2})/6)^{1/2} . Then solve the equation + // 0 = det(\lambda*I - B) = \lambda^{3} - 3*\lambda - det(B) + // which has the solution + // \lambda = 2*cos(1/3 * acos(det(B)/2) +2/3*pi*k ) ; k = 0,1,2 + // when substituting \lambda = 2.cos(theta) and using trig identities. + const Number tr_T = trace(T); + const Number q = tr_T/3.0; + const Number tmp1 = ( T[0][0] - q)*(T[0][0] - q) + + (T[1][1] - q)*(T[1][1] - q) + + (T[2][2] - q)*(T[2][2] - q) + + 2.0 * upp_tri_sq; + const Number p = std::sqrt(tmp1/6.0); + const SymmetricTensor<2,3,Number> B = Number(1.0/p)*(T - q*unit_symmetric_tensor<3,Number>()); + const Number tmp_2 = determinant(B)/2.0; + + // The value of tmp_2 should be within [-1,1], however + // floating point errors might place it slightly outside + // this range. It is therefore necessary to correct for it + const Number phi = + (tmp_2 <= -1.0 ? Number(M_PI/3.0) : + (tmp_2 >= 1.0 ? Number(0.0) : + std::acos(tmp_2)/3.0)); + + // Due to the trigonometric solution, the computed eigenvalues + // should be predictably in the order eig1 >= eig2 >= eig3... + std::array eig_vals + = { { + static_cast(q + 2.0*p*std::cos(phi)), + static_cast(0.0), + static_cast(q + 2.0*p*std::cos(phi + (2.0/3.0*M_PI))) + } + }; + // Use the identity tr(T) = eig1 + eig2 + eig3 + eig_vals[1] = tr_T - eig_vals[0] - eig_vals[2]; + + // ... however, when equal roots exist then floating point + // errors may make this no longer be the case. + // Sort from largest to smallest. + std::sort(eig_vals.begin(), eig_vals.end(), std::greater()); + + return eig_vals; + } +} + + + +namespace internal +{ + namespace SymmetricTensor + { + template + void + tridiagonalize (const dealii::SymmetricTensor<2,dim,Number> &A, + dealii::Tensor<2,dim,Number> &Q, + std::array &d, + std::array &e) + { + // Create some intermediate storage + Number h,g,omega_inv,K,f; + + // Initialize the transformation matrix as the + // identity tensor + Q = dealii::unit_symmetric_tensor(); + + // Make the first row and column to be of the + // desired form + h = 0.0; + for (int i=1; i < dim; i++) + h += A[0][i]*A[0][i]; + + g = 0.0; + if (A[0][1] > 0.0) + g = -std::sqrt(h); + else + g = std::sqrt(h); + e[0] = g; + + std::array u; + for (int i=1; i < dim; i++) + { + u[i] = A[0][i]; + if (i == 1) + u[i] -= g; + } + + std::array q; + const Number omega = h - g * A[0][1]; + if (omega > 0.0) + { + omega_inv = 1.0 / omega; + K = 0.0; + for (int i=1; i < dim; i++) + { + f = 0.0; + for (int j=1; j < dim; j++) + f += A[i][j] * u[j]; + q[i] = omega_inv * f; + K += u[i] * f; + } + K *= 0.5*omega_inv*omega_inv; + + for (int i=1; i < dim; i++) + q[i] = q[i] - K * u[i]; + + d[0] = A[0][0]; + for (int i=1; i < dim; i++) + d[i] = A[i][i] - 2.0*q[i]*u[i]; + + // Store inverse Householder transformation + // in Q + for (int j=1; j < dim; j++) + { + f = omega_inv * u[j]; + for (int i=1; i < dim; i++) + Q[i][j] = Q[i][j] - f*u[i]; + } + + // For dim = 3: Calculate updated A[1][2] and + // store it in e[1] + for (int i=1; i < dim-1; i++) + e[i] = A[i][i+1] - q[i]*u[i+1] - u[i]*q[i+1]; + } + else + { + for (int i=0; i < dim; i++) + d[i] = A[i][i]; + + // For dim = 3: + for (int i=1; i < dim-1; i++) + e[i] = A[i][i+1]; + } + } + + + + template + std::array >,dim> + ql_implicit_shifts (const dealii::SymmetricTensor<2,dim,Number> &A) + { + static_assert(numbers::NumberTraits::is_complex == false, + "This implementation of the QL implicit shift algorithm does " + "not support complex numbers"); + + // Transform A to real tridiagonal form by the Householder method: + // The orthogonal matrix effecting the transformation + // this will ultimately store the eigenvectors + dealii::Tensor<2,dim,Number> Q; + // The diagonal elements of the tridiagonal matrix; + // this will ultimately store the eigenvalues + std::array w; + // The off-diagonal elements of the tridiagonal + std::array ee; + tridiagonalize(A, Q, w, ee); + + // Number of iterations + const unsigned int max_n_it = 30; + + // Transfer the off-diagonal entries to an auxiliary array + // The third element is used only as temporary workspace + std::array e; + for (unsigned int i=0; i >,dim> (); + } + + // Calculate the shift.. + g = (w[l+1] - w[l]) / (e[l] + e[l]); + r = std::sqrt(g*g + 1.0); + // .. and then compute g = d_m - k_s for the + // plane rotation (Press2007a eq 11.4.22) + if (g > 0.0) + g = w[m] - w[l] + e[l]/(g + r); + else + g = w[m] - w[l] + e[l]/(g - r); + + // Perform plane rotation, as is done in the + // standard QL algorithm, followed by Givens + // rotations to recover the tridiagonal form + s = c = 1.0; + p = 0.0; + for (int i=m-1; i >= l; i--) + { + f = s * e[i]; + b = c * e[i]; + + // Branch to recover from underflow + if (std::abs(f) > std::abs(g)) + { + c = g / f; + r = std::sqrt(c*c + 1.0); + e[i+1] = f * r; + c *= (s = 1.0/r); + } + else + { + s = f / g; + r = std::sqrt(s*s + 1.0); + e[i+1] = g * r; + s *= (c = 1.0/r); + } + + g = w[i+1] - p; + r = (w[i] - g)*s + 2.0*c*b; + p = s * r; + w[i+1] = g + p; + g = c*r - b; + + // Form the eigenvectors + for (int k=0; k < dim; k++) + { + t = Q[k][i+1]; + Q[k][i+1] = s*Q[k][i] + c*t; + Q[k][i] = c*Q[k][i] - s*t; + } + } + w[l] -= p; + e[l] = g; + e[m] = 0.0; + } + } + + // Structure the data to be outputted + std::array >,dim> eig_vals_vecs; + for (unsigned int e=0; e + std::array >,dim> + jacobi (dealii::SymmetricTensor<2,dim,Number> A) + { + static_assert(numbers::NumberTraits::is_complex == false, + "This implementation of the Jacobi algorithm does " + "not support complex numbers"); + + // Sums of diagonal resp. off-diagonal elements + Number sd, so; + // sin(phi), cos(phi), tan(phi) and temporary storage + Number s, c, t; + // More temporary storage + Number g, h, z, theta; + // Threshold value + Number thresh; + + // Initialize the transformation matrix as the + // identity tensor + dealii::Tensor<2,dim,Number> Q (dealii::unit_symmetric_tensor()); + + // The diagonal elements of the tridiagonal matrix; + // this will ultimately store the eigenvalues + std::array w; + for (int i=0; i < dim; i++) + w[i] = A[i][i]; + + // Calculate (tr(A))^{2} + sd = trace(A); + sd *= sd; + + // Number of iterations + const unsigned int max_n_it = 50; + for (unsigned int it=0; it <= max_n_it; it++) + { + // Test for convergence + so = 0.0; + for (int p=0; p < dim; p++) + for (int q=p+1; q < dim; q++) + so += std::abs(A[p][q]); + if (so == 0.0) + break; + + // Throw if no convergence is achieved within a + // stipulated number of iterations + if (it == max_n_it) + { + AssertThrow(false, ExcMessage("No convergence in iterative Jacobi eigenvector algorithm.")) + return std::array >,dim> (); + } + + // Compute threshold value which dictates whether or + // not a Jacobi rotation is performed + const unsigned int n_it_skip = 4; + if (it < n_it_skip) + thresh = 0.2 * so / (dim*dim); + else + thresh = 0.0; + + // Perform sweep + for (int p=0; p < dim; p++) + for (int q=p+1; q < dim; q++) + { + g = 100.0 * std::abs(A[p][q]); + + // After a given number of iterations the + // rotation is skipped if the off-diagonal + // element is small + if (it > n_it_skip && + std::abs(w[p]) + g == std::abs(w[p]) && + std::abs(w[q]) + g == std::abs(w[q])) + { + A[p][q] = 0.0; + } + else if (std::abs(A[p][q]) > thresh) + { + // Calculate Jacobi transformation + h = w[q] - w[p]; + + // Compute surrogate for angle theta resulting from + // angle transformation and subsequent smallest solution + // of quadratic equation + if (std::abs(h) + g == std::abs(h)) + { + // Prevent overflow for large theta^2. This computation + // is the algebraic equivalent of t = 1/(2*theta). + t = A[p][q] / h; + } + else + { + theta = 0.5 * h / A[p][q]; + if (theta < 0.0) + t = -1.0 / (std::sqrt(1.0 + theta*theta) - theta); + else + t = 1.0 / (std::sqrt(1.0 + theta*theta) + theta); + } + + // Compute trigonometric functions for rotation + // in such a way as to prevent overflow for + // large theta. + c = 1.0/std::sqrt(1.0 + t*t); + s = t * c; + z = t * A[p][q]; + + // Apply Jacobi transformation... + A[p][q] = 0.0; + w[p] -= z; + w[q] += z; + // ... by executing the various rotations in sequence + for (int r=0; r < p; r++) + { + t = A[r][p]; + A[r][p] = c*t - s*A[r][q]; + A[r][q] = s*t + c*A[r][q]; + } + for (int r=p+1; r < q; r++) + { + t = A[p][r]; + A[p][r] = c*t - s*A[r][q]; + A[r][q] = s*t + c*A[r][q]; + } + for (int r=q+1; r < dim; r++) + { + t = A[p][r]; + A[p][r] = c*t - s*A[q][r]; + A[q][r] = s*t + c*A[q][r]; + } + + // Update the eigenvectors + for (int r=0; r < dim; r++) + { + t = Q[r][p]; + Q[r][p] = c*t - s*Q[r][q]; + Q[r][q] = s*t + c*Q[r][q]; + } + } + } + } + + // Structure the data to be outputted + std::array >,dim> eig_vals_vecs; + for (unsigned int e=0; e + std::array >,2> + hybrid (const dealii::SymmetricTensor<2,2,Number> &A) + { + static_assert(numbers::NumberTraits::is_complex == false, + "This implementation of the 2d Hybrid algorithm does " + "not support complex numbers"); + + const unsigned int dim = 2; + + // Calculate eigenvalues + const std::array w = eigenvalues(A); + + std::array >,dim> eig_vals_vecs; + + Number t, u; // Intermediate storage + t = std::abs(w[0]); + for (unsigned int i=1; i t) + t = u; + } + + if (t < 1.0) + u = t; + else + u = t*t; + + // Estimated maximum roundoff error + const Number error = 256.0 * std::numeric_limits::epsilon() * u*u; + + // Store eigenvalues + eig_vals_vecs[0].first = w[0]; + eig_vals_vecs[1].first = w[1]; + + // Compute eigenvectors + // http://www.math.harvard.edu/archive/21b_fall_04/exhibits/2dmatrices/ + // https://math.stackexchange.com/a/1548616 + if (A[1][0] != 0.0) + { + // First eigenvector + eig_vals_vecs[0].second[0] = w[0] - A[1][1]; + eig_vals_vecs[0].second[1] = A[1][0]; + + // Second eigenvector + eig_vals_vecs[1].second[0] = w[1] - A[1][1]; + eig_vals_vecs[1].second[1] = A[1][0]; + } + else + { + // First eigenvector + eig_vals_vecs[0].second[0] = w[0]; + eig_vals_vecs[0].second[1] = 0.0; + + // Second eigenvector + eig_vals_vecs[1].second[0] = 0.0; + eig_vals_vecs[1].second[1] = w[1]; + } + // Normalize + eig_vals_vecs[0].second /= eig_vals_vecs[0].second.norm(); + eig_vals_vecs[1].second /= eig_vals_vecs[1].second.norm(); + + // If vectors are nearly linearly dependent, or if there might have + // been large cancelations in the calculation of A[i][i] - w[0], fall + // back to QL algorithm + if (eig_vals_vecs[0].second * eig_vals_vecs[1].second > error) + { + return ql_implicit_shifts(A); + } + + return eig_vals_vecs; + } + + + + template + std::array >,3> + hybrid (const dealii::SymmetricTensor<2,3,Number> &A) + { + static_assert(numbers::NumberTraits::is_complex == false, + "This implementation of the 3d Hybrid algorithm does " + "not support complex numbers"); + + const unsigned int dim = 3; + Number norm; // Squared norm or inverse norm of current eigenvector + Number t, u; // Intermediate storage + + // Calculate eigenvalues + const std::array w = eigenvalues(A); + + t = std::abs(w[0]); + for (unsigned int i=1; i t) + t = u; + } + + if (t < 1.0) + u = t; + else + u = t*t; + + // Estimated maximum roundoff error + const Number error = 256.0 * std::numeric_limits::epsilon() * u*u; + + // Initialize the transformation matrix as the + // identity tensor + dealii::Tensor<2,dim,Number> Q; + Q[0][1] = A[0][1]*A[1][2] - A[0][2]*A[1][1]; + Q[1][1] = A[0][2]*A[0][1] - A[1][2]*A[0][0]; + Q[2][1] = A[0][1]*A[0][1]; + + // Calculate first eigenvector by the formula + // v[0] = (A - w[0]).e1 x (A - w[0]).e2 + Q[0][0] = Q[0][1] + A[0][2]*w[0]; + Q[1][0] = Q[1][1] + A[1][2]*w[0]; + Q[2][0] = (A[0][0] - w[0]) * (A[1][1] - w[0]) - Q[2][1]; + norm = Q[0][0]*Q[0][0] + Q[1][0]*Q[1][0] + Q[2][0]*Q[2][0]; + + // If vectors are nearly linearly dependent, or if there might have + // been large cancellations in the calculation of A[i][i] - w[0], fall + // back to QL algorithm + // Note that this simultaneously ensures that multiple eigenvalues do + // not cause problems: If w[0] = w[1], then A - w[0] * I has rank 1, + // i.e. all columns of A - w[0] * I are linearly dependent. + if (norm <= error) + { + return ql_implicit_shifts(A); + } + else // This is the standard branch + { + norm = std::sqrt(1.0 / norm); + for (unsigned j=0; j < dim; j++) + Q[j][0] = Q[j][0] * norm; + } + + // Calculate second eigenvector by the formula + // v[1] = (A - w[1]).e1 x (A - w[1]).e2 + Q[0][1] = Q[0][1] + A[0][2]*w[1]; + Q[1][1] = Q[1][1] + A[1][2]*w[1]; + Q[2][1] = (A[0][0] - w[1]) * (A[1][1] - w[1]) - Q[2][1]; + norm = Q[0][1]*Q[0][1] + Q[1][1]*Q[1][1] + Q[2][1]*Q[2][1]; + if (norm <= error) + { + return ql_implicit_shifts(A); + } + else + { + norm = std::sqrt(1.0 / norm); + for (unsigned int j=0; j < dim; j++) + Q[j][1] = Q[j][1] * norm; + } + + // Calculate third eigenvector according to + // v[2] = v[0] x v[1] + Q[0][2] = Q[1][0]*Q[2][1] - Q[2][0]*Q[1][1]; + Q[1][2] = Q[2][0]*Q[0][1] - Q[0][0]*Q[2][1]; + Q[2][2] = Q[0][0]*Q[1][1] - Q[1][0]*Q[0][1]; + + // Structure the data to be outputted + std::array >,dim> eig_vals_vecs; + for (unsigned int e=0; e + +// Required for instantiation of template functions +#ifdef DEAL_II_WITH_TRILINOS +#include "Sacado.hpp" +#include +#endif + #include +#include DEAL_II_NAMESPACE_OPEN diff --git a/source/base/symmetric_tensor.inst.in b/source/base/symmetric_tensor.inst.in index e0b698fb7f..f9e2186903 100644 --- a/source/base/symmetric_tensor.inst.in +++ b/source/base/symmetric_tensor.inst.in @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// Copyright (C) 2016 by the deal.II authors +// Copyright (C) 2016 - 2017 by the deal.II authors // // This file is part of the deal.II library. // @@ -21,6 +21,44 @@ for (deal_II_dimension : DIMENSIONS; number : REAL_SCALARS) template class SymmetricTensor<4,deal_II_dimension,number>; + + namespace internal + \{ + namespace SymmetricTensor + \{ + template + void + tridiagonalize (const dealii::SymmetricTensor<2,deal_II_dimension,number> &, + dealii::Tensor<2,deal_II_dimension,number> &, + std::array &, + std::array &); + + template + std::array >,deal_II_dimension> + ql_implicit_shifts (const dealii::SymmetricTensor<2,deal_II_dimension,number> &); + + template + std::array >,deal_II_dimension> + jacobi (dealii::SymmetricTensor<2,deal_II_dimension,number>); + +#ifdef DEAL_II_WITH_TRILINOS + template + void + tridiagonalize (const dealii::SymmetricTensor<2,deal_II_dimension,Sacado::Fad::DFad > &, + dealii::Tensor<2,deal_II_dimension,Sacado::Fad::DFad > &, + std::array,deal_II_dimension> &, + std::array,deal_II_dimension-1> &); + + template + std::array, Tensor<1,deal_II_dimension,Sacado::Fad::DFad > >,deal_II_dimension> + ql_implicit_shifts (const dealii::SymmetricTensor<2,deal_II_dimension,Sacado::Fad::DFad > &); + + template + std::array, Tensor<1,deal_II_dimension,Sacado::Fad::DFad > >,deal_II_dimension> + jacobi (dealii::SymmetricTensor<2,deal_II_dimension,Sacado::Fad::DFad >); +#endif + \} + \} } for (deal_II_dimension : DIMENSIONS; number : COMPLEX_SCALARS) @@ -31,3 +69,57 @@ for (deal_II_dimension : DIMENSIONS; number : COMPLEX_SCALARS) template class SymmetricTensor<4,deal_II_dimension,number>; } + + +for (number : REAL_SCALARS) +{ + template + std::array + eigenvalues (const SymmetricTensor<2,1,number> &); + + template + std::array + eigenvalues (const SymmetricTensor<2,2,number> &); + + template + std::array + eigenvalues (const SymmetricTensor<2,3,number> &); + +#ifdef DEAL_II_WITH_TRILINOS + template + std::array,1> + eigenvalues (const SymmetricTensor<2,1,Sacado::Fad::DFad > &); + + template + std::array,2> + eigenvalues (const SymmetricTensor<2,2,Sacado::Fad::DFad > &); + + template + std::array,3> + eigenvalues (const SymmetricTensor<2,3,Sacado::Fad::DFad > &); +#endif + + namespace internal + \{ + namespace SymmetricTensor + \{ + template + std::array >,2> + hybrid (const dealii::SymmetricTensor<2,2,number> &); + + template + std::array >,3> + hybrid (const dealii::SymmetricTensor<2,3,number> &A); + +#ifdef DEAL_II_WITH_TRILINOS + template + std::array, Tensor<1,2,Sacado::Fad::DFad > >,2> + hybrid (const dealii::SymmetricTensor<2,2,Sacado::Fad::DFad > &); + + template + std::array, Tensor<1,3,Sacado::Fad::DFad > >,3> + hybrid (const dealii::SymmetricTensor<2,3,Sacado::Fad::DFad > &A); +#endif + \} + \} +} diff --git a/tests/base/symmetric_tensor_41.cc b/tests/base/symmetric_tensor_41.cc index 820e8321d8..67d4aaae5c 100644 --- a/tests/base/symmetric_tensor_41.cc +++ b/tests/base/symmetric_tensor_41.cc @@ -80,7 +80,7 @@ check_vector (const int index, }; void -test_dim_1 (const enum EigenvectorMethod method, +test_dim_1 (const enum SymmetricTensorEigenvectorMethod method, const double e1, const double tol = 1e-12) { const unsigned int dim = 1; @@ -93,7 +93,7 @@ test_dim_1 (const enum EigenvectorMethod method, } void -test_dim_2 (const enum EigenvectorMethod method, +test_dim_2 (const enum SymmetricTensorEigenvectorMethod method, const double e1, Tensor<1,2> v1, const double e2, const double tol = 1e-12) { @@ -118,7 +118,7 @@ test_dim_2 (const enum EigenvectorMethod method, } void -test_dim_3 (const enum EigenvectorMethod method, +test_dim_3 (const enum SymmetricTensorEigenvectorMethod method, const double e1, Tensor<1,3> v1, const double e2, Tensor<1,3> v2, const double e3, const double tol = 1e-12) @@ -155,7 +155,7 @@ test_dim_3 (const enum EigenvectorMethod method, } -void run_tests(const enum EigenvectorMethod method) +void run_tests(const enum SymmetricTensorEigenvectorMethod method) { // Dim = 1 { @@ -207,7 +207,9 @@ void run_tests(const enum EigenvectorMethod method) // Non-diagonal (large difference) deallog.push("Test 2e"); { - const double tol = (method == dealii::ql_implicit_shifts ? 1e-11 : 1e-12); + const double tol = ( + method == SymmetricTensorEigenvectorMethod::ql_implicit_shifts + ? 1e-11 : 1e-12); test_dim_2(method, 7.2956e8, Tensor<1,2>({3,2}), -5.284e3, tol ); @@ -270,7 +272,11 @@ void run_tests(const enum EigenvectorMethod method) // Non-diagonal (1 large difference) deallog.push("Test 3f"); { - const double tol = (method == dealii::hybrid ? 1e-9 : (method == dealii::ql_implicit_shifts ? 1e-10 : 5e-11)); + const double tol = ( + method == SymmetricTensorEigenvectorMethod::hybrid ? + 1e-9 : + (method == SymmetricTensorEigenvectorMethod::ql_implicit_shifts ? + 1e-10 : 5e-11)); test_dim_3(method, 7.2956e8, Tensor<1,3>({3,2,5}), -4.856e3, Tensor<1,3>({-0.2,3,1}), @@ -281,7 +287,11 @@ void run_tests(const enum EigenvectorMethod method) // Non-diagonal (2 large difference) deallog.push("Test 3g"); { - const double tol = (method == dealii::hybrid ? 1e-8 : (method == dealii::ql_implicit_shifts ? 1e-7 : 2.5e-10)); + const double tol = ( + method == SymmetricTensorEigenvectorMethod::hybrid ? + 1e-8 : + (method == SymmetricTensorEigenvectorMethod::ql_implicit_shifts ? + 1e-7 : 2.5e-10)); test_dim_3(method, 9.274e7, Tensor<1,3>({2,-0.7,1.4}), 2.59343, Tensor<1,3>({0.5,-0.22,-1.42}), @@ -298,14 +308,14 @@ int main() initlog(); deallog.push("Hybrid"); - run_tests(dealii::hybrid); + run_tests(SymmetricTensorEigenvectorMethod::hybrid); deallog.pop(); deallog.push("QL"); - run_tests(dealii::ql_implicit_shifts); + run_tests(SymmetricTensorEigenvectorMethod::ql_implicit_shifts); deallog.pop(); // deallog.push("Jacobi"); - run_tests(dealii::jacobi); + run_tests(SymmetricTensorEigenvectorMethod::jacobi); deallog.pop(); }