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_TOOLS_H
21 : #define LIBMESH_TENSOR_TOOLS_H
22 :
23 : // Local includes
24 : #include "libmesh/libmesh_common.h"
25 : #include "libmesh/compare_types.h"
26 :
27 : #ifdef LIBMESH_HAVE_METAPHYSICL
28 : #include "metaphysicl/dualnumber_decl.h"
29 : #endif
30 :
31 : namespace libMesh
32 : {
33 : // Forward declarations
34 : template <typename T> class TypeVector;
35 : template <typename T> class VectorValue;
36 : template <typename T> class TypeTensor;
37 : template <typename T> class TensorValue;
38 : template <unsigned int N, typename T> class TypeNTensor;
39 :
40 : namespace TensorTools
41 : {
42 : // Any tensor-rank-independent code will need to include
43 : // tensor_tools.h, so we define a product/dot-product here, starting
44 : // with the generic case to apply to scalars.
45 : // Vector specializations will follow.
46 :
47 : template <typename T, typename T2>
48 : inline
49 : typename boostcopy::enable_if_c<ScalarTraits<T>::value && ScalarTraits<T2>::value,
50 : typename CompareTypes<T, T2>::supertype>::type
51 1738360 : inner_product(const T & a, const T2& b)
52 14129248 : { return a * b; }
53 :
54 : template <typename T, typename T2>
55 : inline
56 : typename CompareTypes<T, T2>::supertype
57 110194 : inner_product(const TypeVector<T> & a, const TypeVector<T2> & b)
58 123004 : { return a * b; }
59 :
60 : template <typename T, typename T2>
61 : inline
62 : typename CompareTypes<T, T2>::supertype
63 0 : inner_product(const TypeTensor<T> & a, const TypeTensor<T2> & b)
64 0 : { return a.contract(b); }
65 :
66 : template <unsigned int N, typename T, typename T2>
67 : inline
68 : typename CompareTypes<T, T2>::supertype
69 : inner_product(const TypeNTensor<N,T> & a, const TypeNTensor<N,T2> & b)
70 : { return a.contract(b); }
71 :
72 : template<typename T>
73 : inline
74 2772 : auto norm(const T & a) -> decltype(std::abs(a))
75 2772 : { return std::abs(a); }
76 :
77 : template<typename T>
78 : inline
79 : T norm(std::complex<T> a) { return std::abs(a); }
80 :
81 : template <typename T>
82 : inline
83 : auto norm(const TypeVector<T> & a) -> decltype(TensorTools::norm(T()))
84 : {return std::sqrt(a.norm_sq());}
85 :
86 : template <typename T>
87 : inline
88 0 : auto norm(const VectorValue<T> & a) -> decltype(TensorTools::norm(T()))
89 0 : {return std::sqrt(a.norm_sq());}
90 :
91 : template <typename T>
92 : inline
93 : auto norm(const TypeTensor<T> & a) -> decltype(TensorTools::norm(T()))
94 : {return std::sqrt(a.norm_sq());}
95 :
96 : template <typename T>
97 : inline
98 : auto norm(const TensorValue<T> & a) -> decltype(TensorTools::norm(T()))
99 : {return std::sqrt(a.norm_sq());}
100 :
101 :
102 : template<typename T>
103 : inline
104 313936955 : auto norm_sq(const T & a) -> decltype(std::norm(a))
105 1211492956 : { return std::norm(a); }
106 :
107 : template<typename T>
108 : inline
109 : T norm_sq(std::complex<T> a) { return std::norm(a); }
110 :
111 : template <typename T>
112 : inline
113 : auto norm_sq(const TypeVector<T> & a) -> decltype(std::norm(T()))
114 : {return a.norm_sq();}
115 :
116 : template <typename T>
117 : inline
118 1788809 : auto norm_sq(const VectorValue<T> & a) -> decltype(std::norm(T()))
119 3577618 : {return a.norm_sq();}
120 :
121 : template <typename T>
122 : inline
123 : auto norm_sq(const TypeTensor<T> & a) -> decltype(std::norm(T()))
124 : {return a.norm_sq();}
125 :
126 : template <typename T>
127 : inline
128 : auto norm_sq(const TensorValue<T> & a) -> decltype(std::norm(T()))
129 : {return a.norm_sq();}
130 :
131 :
132 : template<typename T>
133 : inline
134 : bool is_zero(const T & a){ return a.is_zero();}
135 :
136 : // Any tensor-rank-independent code will need to include
137 : // tensor_tools.h, so we define rank-increasing and real-to-number type
138 : // conversion functions here, starting with the generic case to apply
139 : // to scalars.
140 : // Tensor(and higher?) specializations will go in the tensor
141 : // header(s).
142 : template <typename T>
143 : struct IncrementRank
144 : {
145 : typedef VectorValue<T> type;
146 : };
147 :
148 : template <typename T>
149 : struct IncrementRank<VectorValue<T>>
150 : {
151 : typedef TensorValue<T> type;
152 : };
153 :
154 :
155 : template <typename T>
156 : struct IncrementRank<TypeVector<T>>
157 : {
158 : typedef TensorValue<T> type;
159 : };
160 :
161 : template <typename T>
162 : struct IncrementRank<TypeTensor<T>>
163 : {
164 : typedef TypeNTensor<3,T> type;
165 : };
166 :
167 :
168 : template <typename T>
169 : struct IncrementRank<TensorValue<T>>
170 : {
171 : typedef TypeNTensor<3,T> type;
172 : };
173 :
174 : template <unsigned int N, typename T>
175 : struct IncrementRank<TypeNTensor<N,T>>
176 : {
177 : typedef TypeNTensor<N+1,T> type;
178 : };
179 :
180 :
181 : // Also need rank-decreasing case
182 : template <typename T>
183 : struct DecrementRank
184 : {
185 : // The default case is typically an error, but for simpler
186 : // templated code we need it to be compatible with Number
187 : // operations...
188 : typedef T type;
189 : };
190 :
191 : template <typename T>
192 : struct DecrementRank<VectorValue<T>>
193 : {
194 : typedef T type;
195 : };
196 :
197 : template <typename T>
198 : struct DecrementRank<TypeVector<T>>
199 : {
200 : typedef T type;
201 : };
202 :
203 : template <typename T>
204 : struct DecrementRank<TensorValue<T>>
205 : {
206 : typedef VectorValue<T> type;
207 : };
208 :
209 : template <typename T>
210 : struct DecrementRank<TypeTensor<T>>
211 : {
212 : typedef VectorValue<T> type;
213 : };
214 :
215 : template <unsigned int N, typename T>
216 : struct DecrementRank<TypeNTensor<N,T>>
217 : {
218 : typedef TypeNTensor<N-1,T> type;
219 : };
220 :
221 : // Handle the complex-valued case
222 : template <typename T>
223 : struct MakeNumber
224 : {
225 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
226 : typedef std::complex<T> type;
227 : #else
228 : typedef T type;
229 : #endif
230 : };
231 :
232 : template <typename T>
233 : struct MakeNumber<std::complex<T>>
234 : {
235 : // Should this be a compile-time error? we shouldn't need to make numbers out of
236 : // numbers, but then again having the typedef below enables more generic
237 : // programming
238 : typedef std::complex<T> type;
239 : };
240 :
241 :
242 : template <typename T>
243 : struct MakeNumber<TypeVector<T>>
244 : {
245 : typedef TypeVector<typename MakeNumber<T>::type > type;
246 : };
247 :
248 : template <typename T>
249 : struct MakeNumber<VectorValue<T>>
250 : {
251 : typedef VectorValue<typename MakeNumber<T>::type > type;
252 : };
253 :
254 : template <typename T>
255 : struct MakeNumber<TypeTensor<T>>
256 : {
257 : typedef TypeTensor<typename MakeNumber<T>::type> type;
258 : };
259 :
260 : template <typename T>
261 : struct MakeNumber<TensorValue<T>>
262 : {
263 : typedef TypeTensor<typename MakeNumber<T>::type> type;
264 : };
265 :
266 : template <unsigned int N, typename T>
267 : struct MakeNumber<TypeNTensor<N,T>>
268 : {
269 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
270 : typedef TypeNTensor<N,std::complex<T>> type;
271 : #else
272 : typedef TypeNTensor<N,T> type;
273 : #endif
274 : };
275 :
276 : // A utility for determining real-valued (e.g. shape function)
277 : // types from corresponding complex-valued types
278 : template <typename T>
279 : struct MakeReal
280 : {
281 : typedef T type;
282 : };
283 :
284 : template <typename T>
285 : struct MakeReal<std::complex<T>>
286 : {
287 : typedef T type;
288 : };
289 :
290 : template <typename T>
291 : struct MakeReal<TypeVector<T>>
292 : {
293 : typedef TypeVector<typename MakeReal<T>::type > type;
294 : };
295 :
296 : template <typename T>
297 : struct MakeReal<VectorValue<T>>
298 : {
299 : typedef VectorValue<typename MakeReal<T>::type > type;
300 : };
301 :
302 : template <typename T>
303 : struct MakeReal<TypeTensor<T>>
304 : {
305 : typedef TypeTensor<typename MakeReal<T>::type> type;
306 : };
307 :
308 : template <typename T>
309 : struct MakeReal<TensorValue<T>>
310 : {
311 : typedef TypeTensor<typename MakeReal<T>::type> type;
312 : };
313 :
314 : template <unsigned int N, typename T>
315 : struct MakeReal<TypeNTensor<N,T>>
316 : {
317 : typedef TypeNTensor<N,typename MakeReal<T>::type> type;
318 : };
319 :
320 : // Needed for ExactSolution to compile
321 : Number curl_from_grad( const VectorValue<Number> & );
322 :
323 : //! Computes the curl of a vector given the gradient of that vector
324 : VectorValue<Number> curl_from_grad( const TensorValue<Number> & grad );
325 :
326 : /*! Place holder needed for ExactSolution to compile. Will compute the
327 : curl of a tensor given the gradient of that tensor. */
328 : TensorValue<Number> curl_from_grad( const TypeNTensor<3,Number> & grad );
329 :
330 : //! Dummy. Divergence of a scalar not defined, but is needed for ExactSolution to compile
331 : Number div_from_grad( const VectorValue<Number> & grad );
332 :
333 : //! Computes the divergence of a vector given the gradient of that vector
334 : Number div_from_grad( const TensorValue<Number> & grad );
335 :
336 : /*! Place holder needed for ExactSolution to compile. Will compute the
337 : divergence of a tensor given the gradient of that tensor. */
338 : VectorValue<Number> div_from_grad( const TypeNTensor<3,Number> & grad );
339 :
340 : /**
341 : * This helper structure is used to determine whether a template class is one of
342 : * our mathematical structures, like TypeVector, TypeTensor and their descendents
343 : */
344 : template <typename T>
345 : struct MathWrapperTraits
346 : {
347 : static constexpr bool value = false;
348 : };
349 :
350 : template <typename T>
351 : struct MathWrapperTraits<TypeVector<T>>
352 : {
353 : static constexpr bool value = true;
354 : };
355 :
356 : template <typename T>
357 : struct MathWrapperTraits<VectorValue<T>>
358 : {
359 : static constexpr bool value = true;
360 : };
361 :
362 : template <typename T>
363 : struct MathWrapperTraits<TypeTensor<T>>
364 : {
365 : static constexpr bool value = true;
366 : };
367 :
368 : template <typename T>
369 : struct MathWrapperTraits<TensorValue<T>>
370 : {
371 : static constexpr bool value = true;
372 : };
373 :
374 : template <unsigned int N, typename T>
375 : struct MathWrapperTraits<TypeNTensor<N,T>>
376 : {
377 : static constexpr bool value = true;
378 : };
379 :
380 : template <typename T, typename enable = void> struct MakeBaseNumber {};
381 :
382 : template <typename T>
383 : struct MakeBaseNumber<T, typename std::enable_if<ScalarTraits<T>::value>::type> {
384 : typedef typename MakeNumber<T>::type type;
385 : };
386 :
387 : template <template <typename> class Wrapper, typename T>
388 : struct MakeBaseNumber<
389 : Wrapper<T>,
390 : typename std::enable_if<MathWrapperTraits<Wrapper<T>>::value>::type> {
391 : typedef typename MakeBaseNumber<T>::type type;
392 : };
393 :
394 : template <typename T, typename Enable = void>
395 : struct TensorTraits
396 : {
397 : static_assert(always_false<T>,
398 : "Instantiating the generic template of TensorTraits. You must specialize "
399 : "TensorTraits for your type.");
400 : static constexpr unsigned char rank = 0;
401 : };
402 :
403 : template <typename T>
404 : struct TensorTraits<T, typename std::enable_if<ScalarTraits<T>::value>::type>
405 : {
406 : static constexpr unsigned char rank = 0;
407 : };
408 :
409 : template <typename T>
410 : struct TensorTraits<TypeVector<T>>
411 : {
412 : static constexpr unsigned char rank = 1;
413 : };
414 :
415 : template <typename T>
416 : struct TensorTraits<VectorValue<T>>
417 : {
418 : static constexpr unsigned char rank = 1;
419 : };
420 :
421 : template <typename T>
422 : struct TensorTraits<TypeTensor<T>>
423 : {
424 : static constexpr unsigned char rank = 2;
425 : };
426 :
427 : template <typename T>
428 : struct TensorTraits<TensorValue<T>>
429 : {
430 : static constexpr unsigned char rank = 2;
431 : };
432 :
433 : template <typename T, unsigned int N>
434 : struct TensorTraits<TypeNTensor<N, T>>
435 : {
436 : static constexpr unsigned char rank = static_cast<unsigned char>(N);
437 : };
438 :
439 : }//namespace TensorTools
440 :
441 : }//namespace libMesh
442 :
443 : #endif // LIBMESH_TENSOR_TOOLS_H
|
