Line data    Source code

       1              : // Copyright 2024, UChicago Argonne, LLC
       2              : // All Rights Reserved
       3              : // Software Name: NEML2 -- the New Engineering material Model Library, version 2
       4              : // By: Argonne National Laboratory
       5              : // OPEN SOURCE LICENSE (MIT)
       6              : //
       7              : // Permission is hereby granted, free of charge, to any person obtaining a copy
       8              : // of this software and associated documentation files (the "Software"), to deal
       9              : // in the Software without restriction, including without limitation the rights
      10              : // to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
      11              : // copies of the Software, and to permit persons to whom the Software is
      12              : // furnished to do so, subject to the following conditions:
      13              : //
      14              : // The above copyright notice and this permission notice shall be included in
      15              : // all copies or substantial portions of the Software.
      16              : //
      17              : // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
      18              : // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
      19              : // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
      20              : // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
      21              : // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
      22              : // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
      23              : // THE SOFTWARE.
      24              : 
      25              : #include "neml2/tensors/functions/linalg/dsptrf.h"
      26              : #include "neml2/tensors/functions/where.h"
      27              : 
      28              : #include "neml2/tensors/Vec.h"
      29              : #include "neml2/tensors/R2.h"
      30              : #include "neml2/tensors/SR2.h"
      31              : #include "neml2/tensors/SSR4.h"
      32              : #include "neml2/tensors/R4.h"
      33              : 
      34              : namespace neml2::linalg
      35              : {
      36              : SSR4
      37            5 : dsptrf(const Vec & evals, const R2 & evecs, const Vec & f, const Vec & df)
      38              : {
      39              :   // I'm not using `auto` here on purpose to make sure the tensor symmetries are correct.
      40              : 
      41              :   // Helper lambda to handle the degenerate cases
      42           15 :   auto theta = [&](Size i, Size j) -> Scalar
      43              :   {
      44              :     return where(
      45           15 :         evals(i) != evals(j), 0.5 * (f(i) - f(j)) / (evals(i) - evals(j)), 0.25 * (df(i) + df(j)));
      46            5 :   };
      47              : 
      48              :   // Eigenvectors
      49            5 :   const Vec v1 = evecs.col(0);
      50            5 :   const Vec v2 = evecs.col(1);
      51            5 :   const Vec v3 = evecs.col(2);
      52              : 
      53              :   // Some useful tensor products
      54            5 :   const SR2 M_11 = v1.self_outer();
      55            5 :   const SR2 M_22 = v2.self_outer();
      56            5 :   const SR2 M_33 = v3.self_outer();
      57            5 :   const R2 M_12 = v1.outer(v2);
      58            5 :   const R2 M_13 = v1.outer(v3);
      59            5 :   const R2 M_23 = v2.outer(v3);
      60              : 
      61              :   // Derivative contribution from eigenvalues
      62            5 :   const SSR4 P_1 = (df(0) * M_11.outer(M_11) + df(1) * M_22.outer(M_22) + df(2) * M_33.outer(M_33));
      63              : 
      64              :   // Derivative contribution from eigenvectors
      65              :   // Note the symmetrization is handled by the constructor of SSR4(R4)
      66            5 :   const Scalar theta_12 = theta(0, 1);
      67            5 :   const Scalar theta_23 = theta(1, 2);
      68            5 :   const Scalar theta_13 = theta(0, 2);
      69           10 :   const SSR4 P_2 = theta_12 * SSR4(M_12.outer(M_12)) + theta_13 * SSR4(M_13.outer(M_13)) +
      70           15 :                    theta_23 * SSR4(M_23.outer(M_23));
      71              : 
      72           10 :   return P_1 + 4 * P_2;
      73            5 : }
      74              : } // namespace neml2::linalg