Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : #ifndef LIBMESH_TENSOR_VALUE_H
21 : #define LIBMESH_TENSOR_VALUE_H
22 :
23 : // Local includes
24 : #include "libmesh/type_tensor.h"
25 : #include "libmesh/libmesh_device.h"
26 : #include "libmesh/libmesh.h" // for pi
27 :
28 : #ifdef LIBMESH_HAVE_METAPHYSICL
29 : #include "metaphysicl/raw_type.h"
30 : #endif
31 :
32 : namespace libMesh
33 : {
34 :
35 : /**
36 : * This class defines a tensor in LIBMESH_DIM dimensional Real or Complex
37 : * space. The typedef RealTensorValue always defines a real-valued tensor,
38 : * and NumberTensorValue defines a real or complex-valued tensor depending
39 : * on how the library was configured.
40 : *
41 : * \author Roy H. Stogner
42 : * \date 2004
43 : */
44 : template <typename T>
45 1074205258 : class TensorValue : public TypeTensor<T>
46 : {
47 : public:
48 : typedef T value_type;
49 :
50 : template <typename T2>
51 : struct rebind
52 : {
53 : typedef TensorValue<T2> other;
54 : };
55 :
56 : /**
57 : * Empty constructor.
58 : * Gives the tensor 0 in \p LIBMESH_DIM dimensional T space.
59 : */
60 : TensorValue ();
61 :
62 : /**
63 : * Constructor-from-T. By default sets higher dimensional
64 : * entries to 0.
65 : */
66 : explicit TensorValue (const T & xx,
67 : const T & xy=0,
68 : const T & xz=0,
69 : const T & yx=0,
70 : const T & yy=0,
71 : const T & yz=0,
72 : const T & zx=0,
73 : const T & zy=0,
74 : const T & zz=0);
75 :
76 : /**
77 : * Constructor-from-scalars. By default sets higher dimensional
78 : * entries to 0.
79 : */
80 : template <typename Scalar>
81 : explicit TensorValue (const Scalar & xx,
82 1714978888 : const Scalar & xy=0,
83 1714978888 : const Scalar & xz=0,
84 1714978888 : const Scalar & yx=0,
85 1714978888 : const Scalar & yy=0,
86 1714978888 : const Scalar & yz=0,
87 1714978888 : const Scalar & zx=0,
88 1714978888 : const Scalar & zy=0,
89 : typename
90 : std::enable_if<ScalarTraits<Scalar>::value,
91 857489444 : const Scalar>::type & zz=0);
92 :
93 : /**
94 : * Constructor. Takes 1 row vector for LIBMESH_DIM=1
95 : */
96 : template <typename T2>
97 : LIBMESH_DEVICE_INLINE
98 : TensorValue (const TypeVector<T2> & vx);
99 :
100 : /**
101 : * Constructor. Takes 2 row vectors for LIBMESH_DIM=2
102 : */
103 : template <typename T2>
104 : LIBMESH_DEVICE_INLINE
105 : TensorValue (const TypeVector<T2> & vx,
106 : const TypeVector<T2> & vy);
107 :
108 : /**
109 : * Constructor. Takes 3 row vectors for LIBMESH_DIM=3
110 : */
111 : template <typename T2>
112 : LIBMESH_DEVICE_INLINE
113 : TensorValue (const TypeVector<T2> & vx,
114 : const TypeVector<T2> & vy,
115 : const TypeVector<T2> & vz);
116 :
117 :
118 : /**
119 : * Copy-constructor.
120 : */
121 : template <typename T2>
122 : TensorValue (const TensorValue<T2> & p);
123 :
124 : /**
125 : * Copy-constructor.
126 : */
127 : template <typename T2>
128 : TensorValue (const TypeTensor<T2> & p);
129 :
130 :
131 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
132 : /**
133 : * Constructor that takes two \p TypeTensor<Real>
134 : * representing the real and imaginary part as
135 : * arguments.
136 : */
137 : TensorValue (const TypeTensor<Real> & p_re,
138 : const TypeTensor<Real> & p_im);
139 : #endif
140 :
141 : /**
142 : * Assignment-from-scalar operator. Used only to zero out tensors.
143 : */
144 : template <typename Scalar>
145 : LIBMESH_DEVICE_INLINE
146 : typename std::enable_if<
147 : ScalarTraits<Scalar>::value,
148 : TensorValue &>::type
149 418034 : operator = (const Scalar & libmesh_dbg_var(p) )
150 418034 : { libmesh_assert_equal_to (p, Scalar(0)); this->zero(); return *this; }
151 :
152 : /**
153 : * Generate the intrinsic rotation matrix associated with the provided Euler angles. An intrinsic
154 : * rotation leaves bodies in the domain fixed while rotating the coordinate axes. We follow the
155 : * convention described at http://mathworld.wolfram.com/EulerAngles.html (equations 6-14 give the
156 : * entries of the composite transformation matrix). The rotations are performed sequentially about
157 : * the z, x', and z'' axes, in that order. A positive angle for a given step in the rotation
158 : * sequences gives the appearance of rotating an entity in the domain counter-clockwise around the
159 : * rotation axis, although in fact it is the coordinate axes themselves that are rotating. In
160 : * order to give the appearance of a body rotating counter-clockwise, we actually rotate the
161 : * coordinate axes by the \emph negative of the angle passed into the method. All angles should be
162 : * provided in degrees
163 : * @param phi The \emph negative of the angle we will rotate the coordinate axes around the
164 : * original z-axis
165 : * @param theta The \emph negative of the angle we will rotate the coordinate axes around the
166 : * "current" x-axis (post phi), e.g. x'
167 : * @param psi The \emph negative of the angle we will rotate the coordinate axes around the
168 : * "current" z-axis (post phi and theta), e.g. z''
169 : * @return The associated rotation matrix
170 : */
171 : static TensorValue<Real> intrinsic_rotation_matrix(Real phi, Real theta, Real psi);
172 :
173 : /**
174 : * Invert the rotation that would occur if the same angles were provided to \p
175 : * intrinsic_rotation_matrix, e.g. return to the original starting point. All angles should be
176 : * provided in degrees
177 : */
178 : static TensorValue<Real> inverse_intrinsic_rotation_matrix(Real phi, Real theta, Real psi);
179 :
180 : /**
181 : * Generate the extrinsic rotation matrix associated with the provided Euler angles. An extrinsic
182 : * rotation rotates the bodies in the domain and leaves the coordinate axes fixed. We follow the
183 : * convention described at https://en.wikipedia.org/wiki/Euler_angles#Rotation_matrix and we use
184 : * the matrix described by the 'Proper Euler Angles' column and Z1 X2 Z3 row, which indicates that
185 : * the rotations are performed sequentially about the z, x, and z axes, in that order. A
186 : * positive angle yields a counter-clockwise rotation about the axis in question. Note that angles
187 : * should be provided in degrees
188 : * @param angle1_deg rotation angle around z-axis
189 : * @param angle2_deg rotation angle around x-axis (post angle1)
190 : * @param angle3_deg rotation angle around z-axis (post angle1 and angle2)
191 : * @return The associated rotation matrix
192 : */
193 : static TensorValue<Real>
194 : extrinsic_rotation_matrix(Real angle1_deg, Real angle2_deg, Real angle3_deg);
195 :
196 : /**
197 : * Invert the rotation that would occur if the same angles were provided to \p
198 : * extrinsic_rotation_matrix, e.g. return to the original starting point
199 : */
200 : static TensorValue<Real>
201 : inverse_extrinsic_rotation_matrix(Real angle1_deg, Real angle2_deg, Real angle3_deg);
202 : };
203 :
204 : /**
205 : * Useful typedefs to allow transparent switching
206 : * between Real and Complex data types.
207 : */
208 : typedef TensorValue<Real> RealTensorValue;
209 : typedef TensorValue<Number> NumberTensorValue;
210 : typedef RealTensorValue RealTensor;
211 : typedef NumberTensorValue Tensor;
212 :
213 :
214 :
215 : //------------------------------------------------------
216 : // Inline functions
217 : template <typename T>
218 : LIBMESH_DEVICE_INLINE
219 4381144 : TensorValue<T>::TensorValue () :
220 4381144 : TypeTensor<T> ()
221 : {
222 4381144 : }
223 :
224 :
225 :
226 : template <typename T>
227 : LIBMESH_DEVICE_INLINE
228 1513309 : TensorValue<T>::TensorValue (const T & xx,
229 : const T & xy,
230 : const T & xz,
231 : const T & yx,
232 : const T & yy,
233 : const T & yz,
234 : const T & zx,
235 : const T & zy,
236 : const T & zz) :
237 1513309 : TypeTensor<T> (xx,xy,xz,yx,yy,yz,zx,zy,zz)
238 : {
239 1513309 : }
240 :
241 :
242 : template <typename T>
243 : template <typename Scalar>
244 : LIBMESH_DEVICE_INLINE
245 860178844 : TensorValue<T>::TensorValue (const Scalar & xx,
246 : const Scalar & xy,
247 : const Scalar & xz,
248 : const Scalar & yx,
249 : const Scalar & yy,
250 : const Scalar & yz,
251 : const Scalar & zx,
252 : const Scalar & zy,
253 : typename
254 : std::enable_if<ScalarTraits<Scalar>::value,
255 : const Scalar>::type & zz) :
256 860178844 : TypeTensor<T> (xx,xy,xz,yx,yy,yz,zx,zy,zz)
257 : {
258 0 : }
259 :
260 :
261 :
262 : template <typename T>
263 : template <typename T2>
264 : LIBMESH_DEVICE_INLINE
265 : TensorValue<T>::TensorValue (const TensorValue<T2> & p) :
266 : TypeTensor<T> (p)
267 : {
268 : }
269 :
270 :
271 :
272 : template <typename T>
273 : template <typename T2>
274 : LIBMESH_DEVICE_INLINE
275 0 : TensorValue<T>::TensorValue (const TypeVector<T2> & vx) :
276 0 : TypeTensor<T> (vx)
277 : {
278 0 : }
279 :
280 :
281 :
282 : template <typename T>
283 : template <typename T2>
284 : LIBMESH_DEVICE_INLINE
285 : TensorValue<T>::TensorValue (const TypeVector<T2> & vx,
286 : const TypeVector<T2> & vy) :
287 : TypeTensor<T> (vx, vy)
288 : {
289 : }
290 :
291 :
292 :
293 : template <typename T>
294 : template <typename T2>
295 : LIBMESH_DEVICE_INLINE
296 563596 : TensorValue<T>::TensorValue (const TypeVector<T2> & vx,
297 : const TypeVector<T2> & vy,
298 : const TypeVector<T2> & vz) :
299 563596 : TypeTensor<T> (vx, vy, vz)
300 : {
301 563596 : }
302 :
303 :
304 :
305 : template <typename T>
306 : template <typename T2>
307 : LIBMESH_DEVICE_INLINE
308 25101031 : TensorValue<T>::TensorValue (const TypeTensor<T2> & p) :
309 31341757 : TypeTensor<T> (p)
310 : {
311 862155 : }
312 :
313 :
314 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
315 : template <typename T>
316 : LIBMESH_DEVICE_INLINE
317 : TensorValue<T>::TensorValue (const TypeTensor<Real> & p_re,
318 : const TypeTensor<Real> & p_im) :
319 : TypeTensor<T> (Complex (p_re(0,0), p_im(0,0)),
320 : Complex (p_re(0,1), p_im(0,1)),
321 : Complex (p_re(0,2), p_im(0,2)),
322 : Complex (p_re(1,0), p_im(1,0)),
323 : Complex (p_re(1,1), p_im(1,1)),
324 : Complex (p_re(1,2), p_im(1,2)),
325 : Complex (p_re(2,0), p_im(2,0)),
326 : Complex (p_re(2,1), p_im(2,1)),
327 : Complex (p_re(2,2), p_im(2,2)))
328 : {
329 : }
330 : #endif
331 :
332 : template <typename T>
333 : TensorValue<Real>
334 164386 : TensorValue<T>::intrinsic_rotation_matrix(const Real phi, const Real theta, const Real psi)
335 : {
336 : #if LIBMESH_DIM == 3
337 : // We apply a negative sign here or else we don't get the appearance of
338 : // counter-clockwise/right-hand-rule rotation of the bodies with respect to the coordinate axes
339 : // (but as explained in the method doxygen we are *actually* rotating the coordinate axes while
340 : // leaving the bodies fixed)
341 164386 : const Real p = -phi / 180. * pi;
342 164386 : const Real t = -theta / 180. * pi;
343 164386 : const Real s = -psi / 180. * pi;
344 164386 : const Real sp = std::sin(p), cp = std::cos(p);
345 164386 : const Real st = std::sin(t), ct = std::cos(t);
346 164386 : const Real ss = std::sin(s), cs = std::cos(s);
347 :
348 159774 : return TensorValue<Real>(cp * cs - sp * ct * ss,
349 159774 : sp * cs + cp * ct * ss,
350 159774 : st * ss,
351 159774 : -cp * ss - sp * ct * cs,
352 159774 : -sp * ss + cp * ct * cs,
353 159774 : st * cs,
354 159774 : sp * st,
355 164386 : -cp * st,
356 168998 : ct);
357 : #else
358 : libmesh_ignore(phi, theta, psi);
359 : libmesh_error_msg("TensorValue<T>::intrinsic_rotation_matrix() requires libMesh to be compiled "
360 : "with LIBMESH_DIM==3");
361 : // We'll never get here
362 : return TensorValue<Real>();
363 : #endif
364 : }
365 :
366 : template <typename T>
367 : TensorValue<Real>
368 : TensorValue<T>::inverse_intrinsic_rotation_matrix(const Real phi, const Real theta, const Real psi)
369 : {
370 : // The inverse of a rotation matrix is just the transpose
371 : return TensorValue<T>::intrinsic_rotation_matrix(phi, theta, psi).transpose();
372 : }
373 :
374 : template <typename T>
375 : TensorValue<Real>
376 : TensorValue<T>::extrinsic_rotation_matrix(const Real angle1_deg,
377 : const Real angle2_deg,
378 : const Real angle3_deg)
379 : {
380 : #if LIBMESH_DIM == 3
381 : const auto angle1 = angle1_deg / 180. * pi;
382 : const auto angle2 = angle2_deg / 180. * pi;
383 : const auto angle3 = angle3_deg / 180. * pi;
384 : const auto s1 = std::sin(angle1), c1 = std::cos(angle1);
385 : const auto s2 = std::sin(angle2), c2 = std::cos(angle2);
386 : const auto s3 = std::sin(angle3), c3 = std::cos(angle3);
387 :
388 : return TensorValue<Real>(c1 * c3 - c2 * s1 * s3,
389 : -c1 * s3 - c2 * c3 * s1,
390 : s1 * s2,
391 : c3 * s1 + c1 * c2 * s3,
392 : c1 * c2 * c3 - s1 * s3,
393 : -c1 * s2,
394 : s2 * s3,
395 : c3 * s2,
396 : c2);
397 : #else
398 : libmesh_ignore(angle1_deg, angle2_deg, angle3_deg);
399 : libmesh_error_msg("TensorValue<T>::extrinsic_rotation_matrix() requires libMesh to be compiled "
400 : "with LIBMESH_DIM==3");
401 : // We'll never get here
402 : return TensorValue<Real>();
403 : #endif
404 : }
405 :
406 : template <typename T>
407 : TensorValue<Real>
408 : TensorValue<T>::inverse_extrinsic_rotation_matrix(const Real angle1_deg,
409 : const Real angle2_deg,
410 : const Real angle3_deg)
411 : {
412 : // The inverse of a rotation matrix is just the transpose
413 : return TensorValue<T>::extrinsic_rotation_matrix(angle1_deg, angle2_deg, angle3_deg).transpose();
414 : }
415 :
416 : } // namespace libMesh
417 :
418 : #ifdef LIBMESH_HAVE_METAPHYSICL
419 : namespace MetaPhysicL
420 : {
421 : template <typename T>
422 : struct RawType<libMesh::TensorValue<T>>
423 : {
424 : typedef libMesh::TensorValue<typename RawType<T>::value_type> value_type;
425 :
426 : static value_type value (const libMesh::TensorValue<T> & in)
427 : {
428 : value_type ret;
429 : for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
430 : for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
431 : ret(i,j) = raw_value(in(i,j));
432 :
433 : return ret;
434 : }
435 : };
436 :
437 : template <typename T, typename U>
438 : struct ReplaceAlgebraicType<libMesh::TensorValue<T>, U>
439 : {
440 : typedef U type;
441 : };
442 : }
443 : #endif
444 :
445 : #endif // LIBMESH_TENSOR_VALUE_H
|
