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/crystallography/CrystalGeometry.h"
26 :
27 : #include "neml2/base/Factory.h"
28 : #include "neml2/base/TensorName.h"
29 : #include "neml2/tensors/crystallography.h"
30 : #include "neml2/tensors/Scalar.h"
31 : #include "neml2/tensors/Vec.h"
32 : #include "neml2/tensors/R2.h"
33 : #include "neml2/tensors/SR2.h"
34 : #include "neml2/tensors/WR2.h"
35 : #include "neml2/tensors/MillerIndex.h"
36 : #include "neml2/misc/assertions.h"
37 :
38 : namespace neml2::crystallography
39 : {
40 :
41 : register_NEML2_object(CrystalGeometry);
42 :
43 : OptionSet
44 4 : CrystalGeometry::expected_options()
45 : {
46 4 : OptionSet options = Data::expected_options();
47 :
48 4 : options.doc() =
49 4 : "A Data object storing basic crystallographic information for a given crystal system.";
50 :
51 8 : options.set_buffer<TensorName<R2>>("crystal_class");
52 8 : options.set("crystal_class").doc() = "The set of symmetry operations defining the crystal class.";
53 :
54 8 : options.set_buffer<TensorName<Vec>>("lattice_vectors");
55 12 : options.set("lattice_vectors").doc() =
56 4 : "The three lattice vectors defining the crystal translational symmetry";
57 :
58 8 : options.set_buffer<TensorName<MillerIndex>>("slip_directions");
59 8 : options.set("slip_directions").doc() = "A list of Miller indices defining the slip directions";
60 :
61 8 : options.set_buffer<TensorName<MillerIndex>>("slip_planes");
62 4 : options.set("slip_planes").doc() = "A list of Miller indices defining the slip planes";
63 :
64 4 : return options;
65 0 : }
66 :
67 11 : CrystalGeometry::CrystalGeometry(const OptionSet & options)
68 11 : : CrystalGeometry(options, options.get<Factory *>("_factory"))
69 : {
70 11 : }
71 :
72 11 : CrystalGeometry::CrystalGeometry(const OptionSet & options, Factory * factory)
73 : : CrystalGeometry(options,
74 22 : options.get<TensorName<R2>>("crystal_class").resolve(factory),
75 22 : options.get<TensorName<Vec>>("lattice_vectors").resolve(factory),
76 22 : options.get<TensorName<MillerIndex>>("slip_directions").resolve(factory),
77 55 : options.get<TensorName<MillerIndex>>("slip_planes").resolve(factory))
78 : {
79 11 : }
80 :
81 31 : CrystalGeometry::CrystalGeometry(const OptionSet & options,
82 : const R2 & cclass,
83 : const Vec & lattice_vectors,
84 : const MillerIndex & slip_directions,
85 31 : const MillerIndex & slip_planes)
86 : : CrystalGeometry(options,
87 : cclass,
88 : lattice_vectors,
89 31 : setup_schmid_tensors(lattice_vectors, cclass, slip_directions, slip_planes))
90 : {
91 31 : }
92 :
93 31 : CrystalGeometry::CrystalGeometry(const OptionSet & options,
94 : const R2 & cclass,
95 : const Vec & lattice_vectors,
96 31 : std::tuple<Vec, Vec, Scalar, std::vector<Size>> slip_data)
97 : : Data(options),
98 31 : _sym_ops(cclass),
99 31 : _lattice_vectors(declare_buffer<Vec>("lattice_vectors", lattice_vectors)),
100 31 : _reciprocal_lattice_vectors(declare_buffer<Vec>("reciprocal_lattice_vectors",
101 62 : make_reciprocal_lattice(_lattice_vectors))),
102 124 : _slip_directions(declare_buffer<MillerIndex>("slip_directions", "slip_directions")),
103 124 : _slip_planes(declare_buffer<MillerIndex>("slip_planes", "slip_planes")),
104 31 : _cartesian_slip_directions(
105 62 : declare_buffer<Vec>("cartesian_slip_directions", std::get<0>(slip_data))),
106 62 : _cartesian_slip_planes(declare_buffer<Vec>("cartesian_slip_planes", std::get<1>(slip_data))),
107 62 : _burgers(declare_buffer<Scalar>("burgers", std::get<2>(slip_data))),
108 31 : _slip_offsets(std::get<3>(slip_data)),
109 31 : _A(declare_buffer<R2>("schmid_tensors",
110 62 : (_cartesian_slip_directions / _cartesian_slip_directions.norm())
111 62 : .outer(_cartesian_slip_planes / _cartesian_slip_planes.norm()))),
112 93 : _M(declare_buffer<SR2>("symmetric_schmid_tensors", SR2(_A))),
113 124 : _W(declare_buffer<WR2>("skew_symmetric_schmid_tensors", WR2(_A)))
114 : {
115 31 : }
116 :
117 : Vec
118 2 : CrystalGeometry::a1() const
119 : {
120 6 : return _lattice_vectors.batch_index({0});
121 2 : }
122 :
123 : Vec
124 2 : CrystalGeometry::a2() const
125 : {
126 6 : return _lattice_vectors.batch_index({1});
127 2 : }
128 :
129 : Vec
130 2 : CrystalGeometry::a3() const
131 : {
132 6 : return _lattice_vectors.batch_index({2});
133 2 : }
134 :
135 : Vec
136 2 : CrystalGeometry::b1() const
137 : {
138 6 : return _reciprocal_lattice_vectors.batch_index({0});
139 2 : }
140 :
141 : Vec
142 2 : CrystalGeometry::b2() const
143 : {
144 6 : return _reciprocal_lattice_vectors.batch_index({1});
145 2 : }
146 :
147 : Vec
148 2 : CrystalGeometry::b3() const
149 : {
150 6 : return _reciprocal_lattice_vectors.batch_index({2});
151 2 : }
152 :
153 : Size
154 25 : CrystalGeometry::nslip() const
155 : {
156 25 : return _slip_offsets.back();
157 : }
158 :
159 : Size
160 6 : CrystalGeometry::nslip_groups() const
161 : {
162 : // NOLINTNEXTLINE(*-narrowing-conversions)
163 6 : return _slip_offsets.size() - 1;
164 : }
165 :
166 : Size
167 2 : CrystalGeometry::nslip_in_group(Size i) const
168 : {
169 2 : neml_assert_dbg(i < nslip_groups());
170 2 : return _slip_offsets[i + 1] - _slip_offsets[i];
171 : }
172 :
173 : Vec
174 62 : CrystalGeometry::make_reciprocal_lattice(const Vec & lattice_vectors)
175 : {
176 124 : auto a1 = lattice_vectors.batch_index({0});
177 124 : auto a2 = lattice_vectors.batch_index({1});
178 124 : auto a3 = lattice_vectors.batch_index({2});
179 :
180 434 : Vec rl = Vec(at::stack({a2.cross(a3) / a1.dot(a2.cross(a3)),
181 124 : a3.cross(a1) / a2.dot(a3.cross(a1)),
182 186 : a1.cross(a2) / a3.dot(a1.cross(a2))}));
183 :
184 124 : return rl;
185 496 : }
186 :
187 : Vec
188 62 : CrystalGeometry::miller_to_cartesian(const Vec & A, const MillerIndex & d)
189 : {
190 : // Take advantage that a collection of 3 vectors is a R2
191 62 : return R2(ATensor(A)) * d.reduce().to_vec();
192 : }
193 :
194 : std::tuple<Vec, Vec, Scalar, std::vector<Size>>
195 31 : CrystalGeometry::setup_schmid_tensors(const Vec & A,
196 : const R2 & cls,
197 : const MillerIndex & slip_directions,
198 : const MillerIndex & slip_planes)
199 : {
200 : // We need the reciprocol lattice
201 31 : Vec B = make_reciprocal_lattice(A);
202 :
203 : // List of slip directions and planes needs to be consistent
204 31 : if (slip_directions.batch_sizes() != slip_planes.batch_sizes())
205 0 : neml_assert("Input slip directions and planes must have the same batch sizes");
206 :
207 31 : auto bshape = slip_planes.batch_sizes().concrete();
208 31 : auto nbatch = slip_planes.batch_dim();
209 :
210 : // Loop over each slip system
211 31 : std::vector<ATensor> cartesian_slip_directions;
212 31 : std::vector<ATensor> cartesian_slip_planes;
213 31 : std::vector<ATensor> burgers_vectors;
214 62 : std::vector<Size> offsets = {0};
215 :
216 62 : for (Size i = 0; i < bshape[nbatch - 1]; i++)
217 : {
218 : // Get the cartesian slip plane and direction
219 93 : auto cmd = slip_directions.batch_index({indexing::Ellipsis, i});
220 93 : auto cmp = slip_planes.batch_index({indexing::Ellipsis, i});
221 :
222 : // Get the families of symmetry-equivalent planes and directions
223 31 : auto direction_options = unique_bidirectional(cls, miller_to_cartesian(A, cmd));
224 31 : auto plane_options = unique_bidirectional(cls, miller_to_cartesian(B, cmp));
225 :
226 : // Accept the ones that are perpendicular
227 : // We could do this in a vectorized manner, but I don't think it's worth it as
228 : // this code only runs once
229 31 : Size last = offsets.back();
230 217 : for (Size j = 0; j < direction_options.batch_size(-1).concrete(); j++)
231 : {
232 558 : auto di = direction_options.batch_index({indexing::Ellipsis, j});
233 186 : auto dps = plane_options.dot(di);
234 186 : auto inds = at::where(at::isclose(at::abs(dps), at::scalar_tensor(0.0, dps.dtype()))).front();
235 : // We could very easily vectorize this loop, but again whatever
236 558 : for (Size kk = 0; kk < inds.sizes()[0]; kk++)
237 : {
238 744 : Size k = inds.index({kk}).item<Size>();
239 1116 : auto pi = plane_options.batch_index({indexing::Ellipsis, k});
240 372 : cartesian_slip_directions.push_back(di / di.norm());
241 372 : cartesian_slip_planes.push_back(pi / pi.norm());
242 372 : burgers_vectors.push_back(di.norm());
243 372 : last += 1;
244 372 : }
245 186 : }
246 31 : offsets.push_back(last);
247 31 : }
248 :
249 62 : return std::make_tuple(Vec(at::stack(cartesian_slip_directions)),
250 62 : Vec(at::stack(cartesian_slip_planes)),
251 62 : Scalar(at::stack(burgers_vectors)),
252 62 : offsets);
253 1023 : }
254 :
255 : } // namespace neml2
|
