/**
- * Return the nearest orthogonal matrix using a SVD if the determinant is
- * more than a tolerance away from one. The orthogonalization is done by
+ * Return the nearest orthogonal matrix by
* combining the products of the SVD decomposition: $\mathbf U \mathbf{V}^T$,
* where $\mathbf U$ and $\mathbf V$ are computed from the SVD decomposition:
* $\mathbf U \mathbf S \mathbf V^T$,
* effectively replacing $\mathbf S$ with the identity matrix.
* @param tensor The tensor which to find the closest orthogonal
* tensor to.
- * @param tolerance If the $\text{determinant} - 1$ is smaller than
- * this value, it will just return the current tensor.
- * Otherwise it will return the nearest orthogonal tensor.
* @relatesalso Tensor
*/
template <int dim, typename Number>
Tensor<2, dim, Number>
-project_onto_orthogonal_tensors(const Tensor<2, dim, Number> &tensor,
- const double tolerance)
+project_onto_orthogonal_tensors(const Tensor<2, dim, Number> &tensor)
{
- if (std::abs(determinant(tensor) - 1.0) > tolerance)
- {
- Tensor<2, dim, Number> output_tensor;
- FullMatrix<Number> matrix(dim);
- LAPACKFullMatrix<Number> lapack_matrix(dim);
- LAPACKFullMatrix<Number> result(dim);
+ Tensor<2, dim, Number> output_tensor;
+ FullMatrix<Number> matrix(dim);
+ LAPACKFullMatrix<Number> lapack_matrix(dim);
+ LAPACKFullMatrix<Number> result(dim);
- // todo: find or add dealii functionallity to copy in one step.
- matrix.copy_from(tensor);
- lapack_matrix.copy_from(matrix);
+ // todo: find or add dealii functionallity to copy in one step.
+ matrix.copy_from(tensor);
+ lapack_matrix.copy_from(matrix);
- // now compute the svd of the matrices
- lapack_matrix.compute_svd();
+ // now compute the svd of the matrices
+ lapack_matrix.compute_svd();
- // Use the SVD results to orthogonalize: $U V^T$
- lapack_matrix.get_svd_u().mmult(result, lapack_matrix.get_svd_vt());
+ // Use the SVD results to orthogonalize: $U V^T$
+ lapack_matrix.get_svd_u().mmult(result, lapack_matrix.get_svd_vt());
- // todo: find or add dealii functionallity to copy in one step.
- matrix = result;
- matrix.copy_to(output_tensor);
- return output_tensor;
- }
- return tensor;
+ // todo: find or add dealii functionallity to copy in one step.
+ matrix = result;
+ matrix.copy_to(output_tensor);
+ return output_tensor;
}
deallog << "original tensor A = " << A << std::endl;
deallog << "det A = " << determinant(A) << std::endl;
deallog << "A A^T = " << (A * transpose(A)) << std::endl;
- const Tensor<2, dim, number> Q = project_onto_orthogonal_tensors(A, 1.e-10);
+ const Tensor<2, dim, number> Q = project_onto_orthogonal_tensors(A);
deallog << "projected to Q = " << Q << std::endl;
deallog << "det Q = " << determinant(Q) << std::endl;
deallog << "Q Q^T = " << (Q * transpose(Q)) << std::endl;
test(Tensor<2, 3, double>{{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}});
// already orthogonal
test(Tensor<2, 3, double>{{{1, 0, 0}, {0, 0, 1}, {0, 1, 0}}});
+ // not orthogonal, but det = 1
+ test(Tensor<2, 3, double>{{{1, 2, 0}, {0, 1, 0}, {0, 0, 1}}});
// 2D not orthogonal
test(Tensor<2, 2, double>{{{1, 2}, {3, 4}}});
// 2D already orthogonal
DEAL::det Q = -1.00000
DEAL::Q Q^T = 1.00000 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 1.00000
DEAL::
+DEAL:: Tensor<2, 3, double>
+DEAL::original tensor A = 1.00000 2.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 1.00000
+DEAL::det A = 1.00000
+DEAL::A A^T = 5.00000 2.00000 0.00000 2.00000 1.00000 0.00000 0.00000 0.00000 1.00000
+DEAL::projected to Q = 0.707107 0.707107 0.00000 -0.707107 0.707107 0.00000 0.00000 0.00000 1.00000
+DEAL::det Q = 1.00000
+DEAL::Q Q^T = 1.00000 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 1.00000
+DEAL::
DEAL:: Tensor<2, 2, double>
DEAL::original tensor A = 1.00000 2.00000 3.00000 4.00000
DEAL::det A = -2.00000