adress Wolfgangs comments.

diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h
index 81a06ca789368b13a37825d56af5fa452d7a5650..dd4727d2f2fa0b0c3a97061e2a728d4d1d2d3616 100644 (file)
--- a/include/deal.II/base/tensor.h
+++ b/include/deal.II/base/tensor.h
@@ -2618,31 +2618,29 @@ cofactor(const Tensor<2, dim, Number> &t)
 
 
 /**
- * Return the nearest orthogonal matrix using a SVD if the
- * deteriminant is more than a tolerance away from one.
- * The orthogonalization is done by combining the products
- * of the SVD decomposition: $((U*I)*V^T)^T$, where I is the
- * idententy matrix and $U$ and $V$ are computed from the SVD
- * decomposition: $\mathbf U \cdot \mathbf S \cdot \mathbf V^T$
+ * Return the nearest orthogonal matrix using a SVD if the determinant is
+ * more than a tolerance away from one. The orthogonalization is done by
+ * combining the products of the SVD decomposition: $U V^T$, where
+ * $U$ and $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>
 constexpr Tensor<2, dim, Number>
-orthogonalize(const Tensor<2, dim, Number> &tensor, const double tolerance)
+project_onto_orthogonal_tensors(const Tensor<2, dim, Number> &tensor,
+                                const double                  tolerance)
 {
   if (std::abs(determinant(tensor) - 1.0) > tolerance)
     {
-      LAPACKFullMatrix<Number> identity_matrix(dim);
-      for (size_t i = 0; i < dim; i++)
-        {
-          identity_matrix.set(i, i, 1.);
-        }
-
       Tensor<2, dim, Number>   output_tensor;
       FullMatrix<Number>       matrix(dim);
       LAPACKFullMatrix<Number> lapack_matrix(dim);
       LAPACKFullMatrix<Number> result(dim);
-      LAPACKFullMatrix<Number> result2(dim);
 
       // todo: find or add dealii functionallity to copy in one step.
       matrix.copy_from(tensor);
@@ -2651,12 +2649,11 @@ orthogonalize(const Tensor<2, dim, Number> &tensor, const double tolerance)
       // now compute the svd of the matrices
       lapack_matrix.compute_svd();
 
-      // Use the SVD results to orthogonalize: ((U*I)*V^T)^T
-      lapack_matrix.get_svd_u().mmult(result, identity_matrix);
-      result.mmult(result2, (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 = result2;
+      matrix = result;
       matrix.copy_to(output_tensor);
       return output_tensor;
     }

diff --git a/tests/ad_common_tests/tensor_functions_01.h b/tests/ad_common_tests/tensor_functions_01.h
index 85f4fc853a3123b4de4448855419e369dc4e6a82..d0813fa81ad947fd69fb578db37c91e06c198c9f 100644 (file)
--- a/tests/ad_common_tests/tensor_functions_01.h
+++ b/tests/ad_common_tests/tensor_functions_01.h
@@ -55,15 +55,16 @@ test_tensor()
   const Tensor<2, dim, ADNumberType> C5 = A * a;
   const Tensor<2, dim, ADNumberType> C6 = A / a;
 
-  const ADNumberType                 det_A       = determinant(A);
-  const ADNumberType                 tr_A        = trace(A);
-  const Tensor<2, dim, ADNumberType> A_inv       = invert(A);
-  const Tensor<2, dim, ADNumberType> A_T         = transpose(A);
-  const Tensor<2, dim, ADNumberType> A_adj       = adjugate(A);
-  const Tensor<2, dim, ADNumberType> A_cof       = cofactor(A);
-  const Tensor<2, dim, ADNumberType> A_orth      = orthogonalize(A);
-  const ADNumberType                 A_l1_norm   = l1_norm(A, 1e-8);
-  const ADNumberType                 A_linf_norm = linfty_norm(A);
+  const ADNumberType                 det_A = determinant(A);
+  const ADNumberType                 tr_A  = trace(A);
+  const Tensor<2, dim, ADNumberType> A_inv = invert(A);
+  const Tensor<2, dim, ADNumberType> A_T   = transpose(A);
+  const Tensor<2, dim, ADNumberType> A_adj = adjugate(A);
+  const Tensor<2, dim, ADNumberType> A_cof = cofactor(A);
+  const Tensor<2, dim, ADNumberType> A_orth =
+    project_onto_orthogonal_tensors(A);
+  const ADNumberType A_l1_norm   = l1_norm(A, 1e-8);
+  const ADNumberType A_linf_norm = linfty_norm(A);
 
   const ADNumberType A_ddot_B = double_contract<0, 0, 1, 1>(A, B);
   const Tensor<2, dim, ADNumberType> A_dot_B = contract<1, 0>(A, B);

diff --git a/tests/tensors/constexpr_tensor.cc b/tests/tensors/constexpr_tensor.cc
index d1f2f0905e0cff5d72886e94b580854e1b6e6365..11b026535fc2beb6361cb0cb6e5133794260e319 100644 (file)
--- a/tests/tensors/constexpr_tensor.cc
+++ b/tests/tensors/constexpr_tensor.cc
@@ -181,7 +181,7 @@ main()
     DEAL_II_CONSTEXPR const auto dummy_8 = cofactor(a);
     deallog << "Deteriminant before orthogonalization: " << determinant(a)
             << std::endl;
-    const auto dummy_9 = orthogonalize(a, 1e-8);
+    const auto dummy_9 = project_onto_orthogonal_tensors(a, 1e-8);
     deallog << "Deteriminant after  orthogonalization: " << determinant(dummy_9)
             << std::endl;
     Assert(determinant(dummy_9) - 1. < 1e-8, ExcInternalError());
@@ -194,7 +194,7 @@ main()
 
     deallog << "Deteriminant before orthogonalization: "
             << determinant(non_orthogonal) << std::endl;
-    const auto dummy_10 = orthogonalize(non_orthogonal, 1e-8);
+    const auto dummy_10 = project_onto_orthogonal_tensors(non_orthogonal, 1e-8);
     deallog << "Deteriminant after  orthogonalization: "
             << determinant(dummy_10) << std::endl;
     Assert(determinant(dummy_10) - 1. < 1e-8, ExcInternalError());