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/phase_field_fracture/LinearIsotropicStrainEnergyDensity.h"
26 : #include "neml2/misc/errors.h"
27 : #include "neml2/tensors/Vec.h"
28 : #include "neml2/tensors/SR2.h"
29 : #include "neml2/tensors/SSR4.h"
30 : #include "neml2/tensors/Scalar.h"
31 : #include "neml2/tensors/functions/macaulay.h"
32 : #include "neml2/tensors/functions/heaviside.h"
33 : #include "neml2/tensors/functions/linalg/eigh.h"
34 : #include "neml2/tensors/functions/linalg/ieigh.h"
35 : #include "neml2/tensors/functions/linalg/dsptrf.h"
36 : #include "neml2/base/EnumSelection.h"
37 :
38 : namespace neml2
39 : {
40 : register_NEML2_object(LinearIsotropicStrainEnergyDensity);
41 :
42 : OptionSet
43 2 : LinearIsotropicStrainEnergyDensity::expected_options()
44 : {
45 2 : OptionSet options = ElasticityInterface<StrainEnergyDensity, 2>::expected_options();
46 2 : options.doc() =
47 2 : "Calculates elastic strain energy density based on linear elastic isotropic response";
48 4 : options.set<bool>("define_second_derivatives") = true;
49 :
50 8 : EnumSelection type_selection({"NONE", "SPECTRAL", "VOLDEV"}, "NONE");
51 2 : options.set<EnumSelection>("decomposition") = type_selection;
52 4 : options.set("decomposition").doc() =
53 6 : "Strain energy density decomposition types, options are: " + type_selection.candidates_str();
54 :
55 4 : return options;
56 2 : }
57 :
58 10 : LinearIsotropicStrainEnergyDensity::LinearIsotropicStrainEnergyDensity(const OptionSet & options)
59 : : ElasticityInterface<StrainEnergyDensity, 2>(options),
60 10 : _converter(_constant_types, _need_derivs),
61 20 : _decomposition(options.get<EnumSelection>("decomposition").as<DecompositionType>())
62 : {
63 10 : }
64 :
65 : void
66 30 : LinearIsotropicStrainEnergyDensity::set_value(bool out, bool dout_din, bool d2out_din2)
67 : {
68 30 : switch (_decomposition)
69 : {
70 3 : case DecompositionType::NONE:
71 3 : no_decomposition(out, dout_din, d2out_din2);
72 3 : break;
73 12 : case DecompositionType::VOLDEV:
74 12 : voldev_decomposition(out, dout_din, d2out_din2);
75 12 : break;
76 15 : case DecompositionType::SPECTRAL:
77 15 : spectral_decomposition(out, dout_din, d2out_din2);
78 15 : break;
79 0 : default:
80 0 : throw NEMLException("LinearIsotropicStrainEnergyDensity: Unsupported decomposition type.");
81 : }
82 30 : }
83 :
84 : void
85 3 : LinearIsotropicStrainEnergyDensity::no_decomposition(bool out, bool dout_din, bool d2out_din2)
86 : {
87 3 : const auto [K_and_dK, G_and_dG] = _converter.convert(_constants);
88 3 : const auto & [K, dK] = K_and_dK;
89 3 : const auto & [G, dG] = G_and_dG;
90 3 : const auto etr = SR2(_strain).tr();
91 3 : const auto edev = SR2(_strain).dev();
92 3 : const auto I2 = SR2::identity(_strain.options());
93 3 : const auto I4 = SSR4::identity(_strain.options());
94 3 : const auto J = SSR4::identity_dev(_strain.options());
95 :
96 3 : if (out)
97 : {
98 1 : _psie_active = 0.5 * K * etr * etr + G * edev.inner(edev);
99 1 : _psie_inactive = Scalar::zeros_like(_psie_active);
100 : }
101 3 : if (dout_din)
102 : {
103 1 : _psie_active.d(_strain) = K * etr * I2 + 2 * G * edev;
104 : }
105 3 : if (d2out_din2)
106 : {
107 1 : _psie_active.d(_strain, _strain) = K * I4 + 2 * G * J;
108 : }
109 3 : }
110 :
111 : void
112 12 : LinearIsotropicStrainEnergyDensity::voldev_decomposition(bool out, bool dout_din, bool d2out_din2)
113 : {
114 12 : const auto [K_and_dK, G_and_dG] = _converter.convert(_constants);
115 12 : const auto & [K, dK] = K_and_dK;
116 12 : const auto & [G, dG] = G_and_dG;
117 12 : const auto etr = SR2(_strain).tr();
118 12 : const auto edev = SR2(_strain).dev();
119 12 : const auto I2 = SR2::identity(_strain.options());
120 12 : const auto I4 = SSR4::identity(_strain.options());
121 12 : const auto J = SSR4::identity_dev(_strain.options());
122 :
123 : // Decompose based on the trace of strain
124 12 : const auto etr_pos = macaulay(etr);
125 12 : const auto etr_neg = etr - etr_pos;
126 :
127 12 : if (out)
128 : {
129 4 : _psie_active = 0.5 * K * etr_pos * etr_pos + G * edev.inner(edev);
130 4 : _psie_inactive = 0.5 * K * etr_neg * etr_neg;
131 : }
132 12 : if (dout_din)
133 : {
134 4 : _psie_active.d(_strain) = K * etr_pos * I2 + 2 * G * edev;
135 4 : _psie_inactive.d(_strain) = K * etr_neg * I2;
136 : }
137 12 : if (d2out_din2)
138 : {
139 4 : _psie_active.d(_strain, _strain) = K * heaviside(etr) * I4 + 2 * G * J;
140 4 : _psie_inactive.d(_strain, _strain) = K * heaviside(-etr) * I4;
141 : }
142 12 : }
143 :
144 : void
145 15 : LinearIsotropicStrainEnergyDensity::spectral_decomposition(bool out, bool dout_din, bool d2out_din2)
146 : {
147 15 : const auto [lambda_and_dlambda, G_and_dG] = _converter.convert(_constants);
148 15 : const auto & [lambda, dlambda] = lambda_and_dlambda;
149 15 : const auto & [G, dG] = G_and_dG;
150 15 : const auto etr = SR2(_strain).tr();
151 15 : const auto edev = SR2(_strain).dev();
152 15 : const auto I2 = SR2::identity(_strain.options());
153 15 : const auto I4 = SSR4::identity(_strain.options());
154 :
155 : // Decompose based on the eigenvalues of strain
156 15 : const auto etr_pos = macaulay(etr);
157 15 : const auto etr_neg = etr - etr_pos;
158 15 : const auto [evals, evecs] = linalg::eigh(_strain);
159 15 : const auto evals_pos = macaulay(evals);
160 15 : const auto e_pos = linalg::ieigh(evals_pos, evecs);
161 15 : const auto e_neg = _strain - e_pos;
162 :
163 15 : if (out)
164 : {
165 5 : _psie_active = 0.5 * lambda * etr_pos * etr_pos + G * e_pos.inner(e_pos);
166 5 : _psie_inactive = 0.5 * lambda * etr_neg * etr_neg + G * e_neg.inner(e_neg);
167 : }
168 15 : if (dout_din)
169 : {
170 5 : _psie_active.d(_strain) = lambda * etr_pos * I2 + 2 * G * e_pos;
171 5 : _psie_inactive.d(_strain) = lambda * etr_neg * I2 + 2 * G * e_neg;
172 : }
173 15 : if (d2out_din2)
174 : {
175 5 : const auto P4_pos = linalg::dsptrf(evals, evecs, evals_pos, heaviside(evals));
176 5 : const auto P4_neg = SSR4::identity_sym(_strain.options()) - P4_pos;
177 5 : _psie_active.d(_strain, _strain) = lambda * heaviside(etr) * I4 + 2 * G * P4_pos;
178 5 : _psie_inactive.d(_strain, _strain) = lambda * heaviside(-etr) * I4 + 2 * G * P4_neg;
179 5 : }
180 15 : }
181 :
182 : } // namespace neml2
|
