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