Add a new feature for 2nd order symmetric tensor split

doc/news/changes/minor/20230822TaoJin

@@ -0,0 +1,9 @@
+New: Add two tensor functions for the split of a 2nd-order symmetric tensor
+into a positive part and a negative part based on the signs of the eigenvalues
+obtained from the spectrum decomposition. The function positive_negative_split()
+performs the positive-negative split of the 2nd-order symmetric tensor given as 
+the input. The function positive_negative_projectors() not only performs the split,
+but also provides the derivatives (two fourth-order tensors) of the positive/negative 
+part of the tensor with respect to the input tensor.
+<br>
+(Tao Jin, 2023/08/22)

include/deal.II/base/symmetric_tensor.h

@@ -3521,7 +3521,216 @@ outer_product(const SymmetricTensor<2, dim, Number> &t1,
   return tmp;
 }
 
+/**
+ * Perform a spectrum decomposition of a 2nd-order symmetric tensor \a
+ * original_tensor given as the input argument, \f[ \mathrm{original\_tensor} =
+ * \sum_i \lambda_i \, \boldsymbol{n}_i \otimes \boldsymbol{n}_i, \f] where
+ * $\lambda_i$ is the eigenvalue, and $\boldsymbol{n}_i$ is the corresponding
+ * eigenvector. The outputs are the positive part \a positive_part_tensor and
+ * negative part \a negative_part_tensor of the input tensor,
+ * that is,
+ * \f[
+ *   \mathrm{positive\_part\_tensor} = \sum_i <\lambda_i>_+ \boldsymbol{n}_i
+ * \otimes \boldsymbol{n}_i, \quad \mathrm{negative\_part\_tensor} = \sum_i
+ * <\lambda_i>_- \boldsymbol{n}_i \otimes \boldsymbol{n}_i, \f] where
+ * $<\lambda_i>_+ = \mathrm{max}\{ \lambda_i, 0 \}$ and
+ * $<\lambda_i>_- = \mathrm{min}\{ \lambda_i, 0 \}$. Obviously,
+ * \f[
+ *   \mathrm{positive\_part\_tensor}  + \mathrm{negative\_part\_tensor} =
+ * \mathrm{original\_tensor}. \f]
+ *
+ * @param[in] original_tensor The 2nd-order symmetric tensor to be split into
+ * the positive and negative parts
+ * @param[out] positive_part_tensor The positive part of the input tensor in
+ * which the eigenvalues are positive or zero
+ * @param[out] negative_part_tensor The negative part of the input tensor in
+ * which the eigenvalues are negative or zero
+ *
+ * @relatesalso SymmetricTensor
+ */
+template <int dim, typename Number>
+void
+positive_negative_split(const SymmetricTensor<2, dim, Number> &original_tensor,
+                        SymmetricTensor<2, dim, Number> &positive_part_tensor,
+                        SymmetricTensor<2, dim, Number> &negative_part_tensor)
+{
+  Assert(dim <= 3, ExcMessage("dim should not be larger than 3."));
+
+  std::array<std::pair<Number, Tensor<1, dim, Number>>, dim> eigen_system;
+  std::vector<Number>                                        eigen_values(dim);
+  std::vector<Tensor<1, dim, Number>>                        eigen_vectors(dim);
+
+  eigen_system = eigenvectors(original_tensor);
+
+  for (int i = 0; i < dim; i++)
+    {
+      eigen_values[i]  = eigen_system[i].first;
+      eigen_vectors[i] = eigen_system[i].second;
+    }
+
+  positive_part_tensor = 0;
+  for (int i = 0; i < dim; i++)
+    positive_part_tensor +=
+      std::fmax(eigen_values[i], 0.0) *
+      symmetrize(outer_product(eigen_vectors[i], eigen_vectors[i]));
+
+  negative_part_tensor = 0;
+  for (int i = 0; i < dim; i++)
+    negative_part_tensor +=
+      std::fmin(eigen_values[i], 0.0) *
+      symmetrize(outer_product(eigen_vectors[i], eigen_vectors[i]));
+}
+
+/**
+ * This function is similar to the function positive_negative_split(). That is,
+ * perform a spectrum decomposition of a 2nd-order symmetric tensor \a
+ * original_tensor given as the input argument, and split it into a positive
+ * part \a positive_part_tensor and a negative part \a negative_part_tensor.
+ * Moreover, this function also provides the derivatives. Let $\mathbf{A}$
+ * represent the input 2nd-order symmetric tensor \a original_tensor,
+ * $\mathbf{A}^+$ represent the positive part \a positive_part_tensor, and
+ * $\mathbf{A}^-$ represent the negative part \a negative_part_tensor. Then, two
+ * fourth-order tensors are defined as
+ * \f[
+ *   \mathbb{P}^+ = \frac{\partial \mathbf{A}^+}{\partial \mathbf{A}}, \quad
+ *   \mathbb{P}^- = \frac{\partial \mathbf{A}^-}{\partial \mathbf{A}},
+ * \f]
+ * where $\mathbb{P}^+$ is the \a positive_projector and $\mathbb{P}^-$ is the
+ * \a negative_projector. These two fourth-order tensors satisfy the following
+ * properties: \f[ \mathbb{P}^+ : \mathbf{A} = \mathbf{A}^+, \quad \mathbb{P}^-
+ * : \mathbf{A} = \mathbf{A}^-. \f] Since $\mathbb{P}^+$ and $\mathbb{P}^-$ are
+ * 4th-order projectors, \f[ \mathbb{P}^+ : \mathbf{A}^+ = \mathbf{A}^+, \quad
+ * \mathbb{P}^- : \mathbf{A}^- = \mathbf{A}^-, \quad \mathbb{P}^+ : \mathbf{A}^-
+ * = \mathbb{P}^- : \mathbf{A}^+ = \mathbf{0}. \f] Lastly, \f[ \mathbb{P}^+ +
+ * \mathbb{P}^- = \mathbb{S}, \f] where $\mathbb{S}$ is the fourth-order
+ * symmetric identity tensor Physics::Elasticity::StandardTensors< dim >::S.
+ *
+ * @param[in] original_tensor The 2nd-order symmetric tensor to be split into
+ * the positive and negative parts
+ * @param[out] positive_part_tensor The positive part of the input tensor in
+ * which the eigenvalues are positive or zero
+ * @param[out] negative_part_tensor The negative part of the input tensor in
+ * which the eigenvalues are negative or zero
+ * @param[out] positive_projector The fourth-order positive projection tensor
+ * $\mathbb{P}^+$
+ * @param[out] negative_projector The fourth-order negative projection tensor
+ * $\mathbb{P}^-$
+ *
+ * @relatesalso SymmetricTensor
+ */
+template <int dim, typename Number>
+void
+positive_negative_projectors(
+  const SymmetricTensor<2, dim, Number> &original_tensor,
+  SymmetricTensor<2, dim, Number> &      positive_part_tensor,
+  SymmetricTensor<2, dim, Number> &      negative_part_tensor,
+  SymmetricTensor<4, dim, Number> &      positive_projector,
+  SymmetricTensor<4, dim, Number> &      negative_projector)
+{
+  Assert(dim <= 3, ExcMessage("dim should not be larger than 3."));
+
+  auto heaviside_function{[](const double x) {
+    if (std::fabs(x) < 1.0e-16)
+      return 0.5;
+    if (x > 0)
+      return 1.0;
+    else
+      return 0.0;
+  }};
+
+  std::array<std::pair<Number, Tensor<1, dim, Number>>, dim> eigen_system;
+  std::vector<Number>                                        eigen_values(dim);
+  std::vector<Tensor<1, dim, Number>>                        eigen_vectors(dim);
+
+  eigen_system = eigenvectors(original_tensor);
+
+  for (int i = 0; i < dim; i++)
+    {
+      eigen_values[i]  = eigen_system[i].first;
+      eigen_vectors[i] = eigen_system[i].second;
+    }
 
+  positive_part_tensor = 0;
+  for (int i = 0; i < dim; i++)
+    positive_part_tensor +=
+      std::fmax(eigen_values[i], 0.0) *
+      symmetrize(outer_product(eigen_vectors[i], eigen_vectors[i]));
+
+  negative_part_tensor = 0;
+  for (int i = 0; i < dim; i++)
+    negative_part_tensor +=
+      std::fmin(eigen_values[i], 0.0) *
+      symmetrize(outer_product(eigen_vectors[i], eigen_vectors[i]));
+
+  std::vector<SymmetricTensor<2, dim, Number>> M(dim);
+  for (int a = 0; a < dim; a++)
+    M[a] = symmetrize(outer_product(eigen_vectors[a], eigen_vectors[a]));
+
+  std::vector<SymmetricTensor<4, dim, Number>> Q(dim);
+  for (int a = 0; a < dim; a++)
+    Q[a] = outer_product(M[a], M[a]);
+
+  std::vector<std::vector<SymmetricTensor<4, dim, Number>>> G(
+    dim, std::vector<SymmetricTensor<4, dim, Number>>(dim));
+  for (int a = 0; a < dim; a++)
+    for (int b = 0; b < dim; b++)
+      for (int i = 0; i < dim; i++)
+        for (int j = 0; j < dim; j++)
+          for (int k = 0; k < dim; k++)
+            for (int l = 0; l < dim; l++)
+              G[a][b][i][j][k][l] =
+                M[a][i][k] * M[b][j][l] + M[a][i][l] * M[b][j][k];
+
+  // positive P
+  positive_projector = 0;
+  for (int a = 0; a < dim; a++)
+    {
+      double lambda_a = eigen_values[a];
+      positive_projector += heaviside_function(lambda_a) * Q[a];
+      for (int b = 0; b < dim; b++)
+        {
+          if (b != a)
+            {
+              double lambda_b = eigen_values[b];
+
+              double v_ab = 0.0;
+              if (std::fabs(lambda_a - lambda_b) > 1.0e-12)
+                v_ab = (std::fmax(lambda_a, 0.0) - std::fmax(lambda_b, 0.0)) /
+                       (lambda_a - lambda_b);
+              else
+                v_ab = 0.5 * (heaviside_function(lambda_a) +
+                              heaviside_function(lambda_b));
+
+              positive_projector += 0.5 * v_ab * 0.5 * (G[a][b] + G[b][a]);
+            }
+        }
+    }
+
+  // negative P
+  negative_projector = 0;
+  for (int a = 0; a < dim; a++)
+    {
+      double lambda_a = eigen_values[a];
+      negative_projector += heaviside_function(-lambda_a) * Q[a];
+      for (int b = 0; b < dim; b++)
+        {
+          if (b != a)
+            {
+              double lambda_b = eigen_values[b];
+
+              double v_ab = 0.0;
+              if (std::fabs(lambda_a - lambda_b) > 1.0e-12)
+                v_ab = (std::fmin(lambda_a, 0.0) - std::fmin(lambda_b, 0.0)) /
+                       (lambda_a - lambda_b);
+              else
+                v_ab = 0.5 * (heaviside_function(-lambda_a) +
+                              heaviside_function(-lambda_b));
+
+              negative_projector += 0.5 * v_ab * 0.5 * (G[a][b] + G[b][a]);
+            }
+        }
+    }
+}
 
 /**
  * Return the symmetrized version of a full rank-2 tensor, i.e.

tests/physics/positive_negative_split.cc

@@ -0,0 +1,138 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 - 2023 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+/*
+ * Author: Tao Jin
+ *         University of Ottawa, Ottawa, Ontario, Canada
+ *         August 2023
+ *
+ * Test the positive-negative split of a 2nd-order symmetric tensor
+ * and the positive and negative 4th-order projectors
+ */
+
+
+#include <deal.II/base/symmetric_tensor.h>
+
+#include <deal.II/physics/elasticity/standard_tensors.h>
+
+#include "../tests.h"
+
+template <int dim>
+void
+positive_negative_split_test()
+{
+  using namespace dealii;
+
+  SymmetricTensor<2, dim> random_tensor;
+  srand(time(0));
+
+  for (unsigned int i = 0; i < dim; i++)
+    for (unsigned int j = 0; j <= i; j++)
+      {
+        random_tensor[i][j] = ((double)rand() / (RAND_MAX));
+        if (j != i)
+          random_tensor[j][i] = random_tensor[i][j];
+      }
+  SymmetricTensor<2, dim> positive_part_tensor, negative_part_tensor;
+
+  positive_negative_split(random_tensor,
+                          positive_part_tensor,
+                          negative_part_tensor);
+
+  bool positive_negative_split_success = true;
+
+  // test: (A^+) + (A^-) = A
+  if ((positive_part_tensor + negative_part_tensor - random_tensor).norm() >
+      1.0e-12 * random_tensor.norm())
+    positive_negative_split_success = false;
+
+  if (!positive_negative_split_success)
+    Assert(false, ExcMessage("Positive-negative split failed!"));
+
+  SymmetricTensor<4, dim> positive_projector, negative_projector;
+  positive_negative_projectors(random_tensor,
+                               positive_part_tensor,
+                               negative_part_tensor,
+                               positive_projector,
+                               negative_projector);
+
+  bool                    positive_projector_success = true;
+  SymmetricTensor<2, dim> projected_positive_tensor;
+  projected_positive_tensor = positive_projector * random_tensor;
+
+  // test: (P^+) : A = (A^+)
+  if ((projected_positive_tensor - positive_part_tensor).norm() >
+      1.0e-12 * random_tensor.norm())
+    positive_projector_success = false;
+
+  // test: (P^+) : (A^+) = (A^+)
+  if ((positive_projector * projected_positive_tensor - positive_part_tensor)
+        .norm() > 1.0e-12 * random_tensor.norm())
+    positive_projector_success = false;
+
+  // test: (P^+) : (A^-) = 0
+  if ((positive_projector * negative_part_tensor).norm() >
+      1.0e-12 * random_tensor.norm())
+    positive_projector_success = false;
+
+  bool                    negative_projector_success = true;
+  SymmetricTensor<2, dim> projected_negative_tensor;
+  projected_negative_tensor = negative_projector * random_tensor;
+
+  // test: (P^-) : A = (A^-)
+  if ((projected_negative_tensor - negative_part_tensor).norm() >
+      1.0e-12 * random_tensor.norm())
+    negative_projector_success = false;
+
+  // test: (P^-) : (A^-) = (A^-)
+  if ((negative_projector * projected_negative_tensor - negative_part_tensor)
+        .norm() > 1.0e-12 * random_tensor.norm())
+    negative_projector_success = false;
+
+  // test: (P^-) : (A^+) = 0
+  if ((negative_projector * positive_part_tensor).norm() >
+      1.0e-12 * random_tensor.norm())
+    negative_projector_success = false;
+
+  // test: (P^+) + (P^-) = S (S is 4th-order symmetric identity tensor)
+  if ((positive_projector + negative_projector -
+       Physics::Elasticity::StandardTensors<dim>::S)
+        .norm() > 1.0e-12)
+    {
+      positive_projector_success = false;
+      negative_projector_success = false;
+    }
+
+  if (!positive_projector_success)
+    Assert(false, ExcMessage("Positive projector failed!"));
+
+  if (!negative_projector_success)
+    Assert(false, ExcMessage("Negative projector failed!"));
+}
+
+
+int
+main()
+{
+  initlog();
+
+  positive_negative_split_test<1>();
+  positive_negative_split_test<2>();
+  positive_negative_split_test<3>();
+
+  deallog << "OK" << std::endl;
+
+  return 0;
+}

tests/physics/positive_negative_split.output

@@ -0,0 +1,2 @@
+
+DEAL::OK