From bd7d41880ea0bb9f2e9e15d36969b1881ed894ed Mon Sep 17 00:00:00 2001 From: Tao Jin Date: Tue, 22 Aug 2023 18:31:51 -0400 Subject: [PATCH] new feature to split a 2nd-order symmetric tensor into a positive part and a negative part split a 2nd-order symmetric tensor based on the signs of the eigenvalues obtained from a spectrum decomposition, also provide the corresponding 4th-order tensors that are the derivatives of the positive/negative part of the tensor with respect to the original tensor. --- include/deal.II/base/symmetric_tensor.h | 209 +++++++++++++++++++ tests/physics/positive_negative_split.cc | 138 ++++++++++++ tests/physics/positive_negative_split.output | 2 + 3 files changed, 349 insertions(+) create mode 100644 tests/physics/positive_negative_split.cc create mode 100644 tests/physics/positive_negative_split.output diff --git a/include/deal.II/base/symmetric_tensor.h b/include/deal.II/base/symmetric_tensor.h index 646265e8f6..d221be3422 100644 --- a/include/deal.II/base/symmetric_tensor.h +++ b/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 +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>, dim> eigen_system; + std::vector eigen_values(dim); + std::vector> 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 +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>, dim> eigen_system; + std::vector eigen_values(dim); + std::vector> 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> M(dim); + for (int a = 0; a < dim; a++) + M[a] = symmetrize(outer_product(eigen_vectors[a], eigen_vectors[a])); + + std::vector> Q(dim); + for (int a = 0; a < dim; a++) + Q[a] = outer_product(M[a], M[a]); + + std::vector>> G( + dim, std::vector>(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. diff --git a/tests/physics/positive_negative_split.cc b/tests/physics/positive_negative_split.cc new file mode 100644 index 0000000000..e17f1787da --- /dev/null +++ b/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 + +#include + +#include "../tests.h" + +template +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::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; +} diff --git a/tests/physics/positive_negative_split.output b/tests/physics/positive_negative_split.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/physics/positive_negative_split.output @@ -0,0 +1,2 @@ + +DEAL::OK -- 2.39.5