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/models/FischerBurmeister.h"
      26              : #include "neml2/tensors/functions/sqrt.h"
      27              : #include "neml2/tensors/assertions.h"
      28              : 
      29              : namespace neml2
      30              : {
      31              : register_NEML2_object(FischerBurmeister);
      32              : OptionSet
      33            2 : FischerBurmeister::expected_options()
      34              : {
      35            2 :   OptionSet options = Model::expected_options();
      36            2 :   options.doc() = "By default,  if \\f$ a \\ge 0, "
      37              :                   "b \\ge 0, ab = 0 \\f$ then  the Fischer Burmeister (FB) condition is:\\f$"
      38              :                   "a+b-\\sqrt(a^2+b^2)\\f$, where a, b is the first_var and second_var "
      39              :                   "respectively and first_inequality = second_inequality = 'GE'. One could set "
      40              :                   "first_inequality = 'LE' (i.e. \\f$ a \\le 0, "
      41              :                   "b \\ge 0, ab = 0 \\f$, FB conditions is \\f$"
      42            2 :                   "-a+b-\\sqrt(a^2+b^2) \\f$). Same goes for second_inequality = 'LE'.";
      43              : 
      44            6 :   options.set_input("first_var") = VariableName(STATE, "a");
      45            2 :   options.set("first_var").doc() = "First condition";
      46              : 
      47            6 :   options.set_input("second_var") = VariableName(STATE, "b");
      48            4 :   options.set("second_var").doc() = "Second condition";
      49              : 
      50            8 :   EnumSelection conda({"GE", "LE"}, "GE");
      51            2 :   options.set<EnumSelection>("first_inequality") = conda;
      52            4 :   options.set("first_inequality").doc() = "Type of inequality for the first variable."
      53            4 :                                           "Default: GE. Options are " +
      54            6 :                                           conda.candidates_str();
      55              : 
      56            8 :   EnumSelection condb({"GE", "LE"}, "GE");
      57            2 :   options.set<EnumSelection>("second_inequality") = condb;
      58            4 :   options.set("second_inequality").doc() = "Type of inequality for the second variable."
      59            4 :                                            "Default: GE. Options are " +
      60            6 :                                            condb.candidates_str();
      61              : 
      62            6 :   options.set_output("fischer_burmeister") = VariableName(STATE, "fb");
      63            2 :   options.set("fischer_burmeister").doc() = "Fischer Burmeister condition";
      64              : 
      65            4 :   return options;
      66            2 : }
      67              : 
      68            3 : FischerBurmeister::FischerBurmeister(const OptionSet & options)
      69              :   : Model(options),
      70            3 :     _a(declare_input_variable<Scalar>("first_var")),
      71            3 :     _b(declare_input_variable<Scalar>("second_var")),
      72            6 :     _conda(options.get<EnumSelection>("first_inequality")),
      73            3 :     _condb(options.get<EnumSelection>("second_inequality")),
      74            6 :     _fb(declare_output_variable<Scalar>("fischer_burmeister"))
      75              : {
      76            3 : }
      77              : 
      78              : void
      79            6 : FischerBurmeister::set_value(bool out, bool dout_din, bool d2out_din2)
      80              : {
      81            6 :   neml_assert_dbg(!d2out_din2, "Second derivative not implemented.");
      82            6 :   neml_assert_dbg(_a.scalar_type() == _b.scalar_type(),
      83              :                   "First and second variables must have the same scalar type.");
      84              : 
      85            6 :   auto ia = 1.0;
      86           12 :   if (_conda == "LE")
      87            4 :     ia = -1.0;
      88              : 
      89            6 :   auto ib = 1.0;
      90           12 :   if (_condb == "LE")
      91            2 :     ib = -1.0;
      92              : 
      93            6 :   if (out)
      94              :   {
      95            3 :     _fb = _a * ia + _b * ib - sqrt(_a * _a + _b * _b);
      96              :   }
      97              : 
      98            6 :   if (dout_din)
      99              :   {
     100            3 :     const auto eps = machine_precision(_a.scalar_type());
     101            3 :     _fb.d(_a) = ia - _a / sqrt(_a * _a + _b * _b + eps);
     102            3 :     _fb.d(_b) = ib - _b / sqrt(_a * _a + _b * _b + eps);
     103            3 :   }
     104            6 : }
     105              : }