BELFEM 0.9.0
Berkeley Lab Finite Element Framework
Loading...
Searching...
No Matches
cl_Tensor.hpp
Go to the documentation of this file.
1/*
2 * BELFEM -- The Berkeley Lab Finite Element Framework
3 * Copyright (c) 2026, The Regents of the University of California,
4 * through Lawrence Berkeley National Laboratory (subject to receipt of any required
5 * approvals from the U.S. Dept. of Energy). All rights reserved.
6 *
7 * Developers: Christian Messe, Gregory Giard
8 *
9 * See the top-level LICENSE file for the complete license and disclaimer.
10 */
11
12#ifndef BELFEM_CL_TENSOR_HPP
13#define BELFEM_CL_TENSOR_HPP
14
15#include <cstring> // for std::memcpy
16#include <algorithm> // for std::sort std::fill_n
17#include <functional> // for std::plus std::minus
18
19#include "typedefs.hpp"
20#include "assert.hpp"
21#include "cl_Matrix.hpp"
22#include "fn_TR_mat_to_ten.hpp"
23#include "fn_TR_ten_to_mat.hpp"
24#include "fn_TR_contract42.hpp"
25#include "fn_TR_contract44.hpp"
26#include "fn_TR_equal_equal.hpp"
27
28#include "fn_TR_fill.hpp"
29#include "fn_compliance_matrix.hpp"
30
31#include "fn_inv.hpp"
32namespace belfem
33{
41 template< typename T >
42 class Tensor
43 {
44//----------------------------------------------------------------------------
45
46 // indices
47 const index_t mSizeI;
48 const index_t mSizeJ;
49 const index_t mSizeK;
50 const index_t mSizeL;
51
52 // memory jumps
53 const index_t mOffsetJ;
54 const index_t mOffsetK;
55 const index_t mOffsetL;
56
57 const index_t mOrder;
58 const index_t mCapacity;
59
60 T * mData = nullptr;
61
62//----------------------------------------------------------------------------
63 public:
64//----------------------------------------------------------------------------
65
66 // create an empty tensor
67 Tensor( const index_t aSizeI,
68 const index_t aSizeJ,
69 const index_t aSizeK,
70 const index_t aSizeL ) :
71 mSizeI( aSizeI ),
72 mSizeJ( aSizeJ ),
73 mSizeK( aSizeK ),
74 mSizeL( aSizeL ),
75 mOffsetJ( aSizeI ),
76 mOffsetK( aSizeI * aSizeJ ),
77 mOffsetL( aSizeI * aSizeJ * aSizeK ),
78 mOrder( 4 ),
79 mCapacity( aSizeI * aSizeJ * aSizeK * aSizeL )
80 {
81 BELFEM_ASSERT( mCapacity > 0, "Tensor sizes must not be zero" );
82
83 mData = ( T * ) malloc( ( mCapacity ) * sizeof( T ) );
84 }
85
86 // create an empty tensor
87 Tensor( const index_t aSizeI,
88 const index_t aSizeJ,
89 const index_t aSizeK ) :
90 mSizeI( aSizeI ),
91 mSizeJ( aSizeJ ),
92 mSizeK( aSizeK ),
93 mSizeL( 1 ),
94 mOffsetJ( aSizeI ),
95 mOffsetK( aSizeI * aSizeJ ),
96 mOffsetL( aSizeI * aSizeJ * aSizeK ),
97 mOrder( 3 ),
98 mCapacity( aSizeI * aSizeJ * aSizeK )
99 {
100 BELFEM_ASSERT( mCapacity > 0, "Tensor sizes must not be zero" );
101
102 mData = ( T * ) malloc( ( mCapacity ) * sizeof( T ) );
103 }
104
105 // create an empty tensor and initialize it with values
106 Tensor( const index_t aSizeI,
107 const index_t aSizeJ,
108 const index_t aSizeK,
109 const index_t aSizeL,
110 const real aValue ) :
111 Tensor( aSizeI, aSizeJ, aSizeK, aSizeL )
112 {
113 this->fill( aValue );
114 }
115
116 // create a tensor form an elasticity matrix
117 Tensor( const Matrix< real > & aElasticityMatrix ) :
118 Tensor( 3, 3, 3, 3 )
119 {
120 BELFEM_ASSERT( aElasticityMatrix.n_rows() == 6
121 && aElasticityMatrix.n_cols() == 6,
122 "Matrix must be allocated as 6x6" );
123
124 tensor::mat_to_ten( aElasticityMatrix, mData );
125 }
126
127//----------------------------------------------------------------------------
128
129 // copy constructor
130 Tensor( const Tensor< T > & aOther ) :
131 mSizeI( aOther.size_i() ),
132 mSizeJ( aOther.size_j() ),
133 mSizeK( aOther.size_k() ),
134 mSizeL( aOther.size_l() ),
135 mOffsetJ( aOther.mOffsetJ ),
136 mOffsetK( aOther.mOffsetK ),
137 mOffsetL( aOther.mOffsetL ),
138 mOrder( aOther.mOrder ),
139 mCapacity( aOther.mCapacity )
140 {
141 mData = ( T * ) malloc( mCapacity * sizeof( T ) );
142 std::memcpy( mData, aOther.mData, mCapacity * sizeof( T ) );
143 }
144
145//----------------------------------------------------------------------------
146
147 // move constructor
148 Tensor( Tensor< T > && aOther ) :
149 mSizeI( aOther.size_i() ),
150 mSizeJ( aOther.size_j() ),
151 mSizeK( aOther.size_k() ),
152 mSizeL( aOther.size_l() ),
153 mOffsetJ( aOther.mOffsetJ ),
154 mOffsetK( aOther.mOffsetK ),
155 mOffsetL( aOther.mOffsetL ),
156 mOrder( aOther.mOrder ),
157 mCapacity( aOther.mCapacity )
158 {
159 mData = aOther.mData ;
160 aOther.mData = nullptr ;
161 }
162
163//----------------------------------------------------------------------------
164
165 ~Tensor()
166 {
167 free( mData );
168 }
169
170//------------------------------------------------------------------------------
171// SIZES
172//------------------------------------------------------------------------------
173
174 index_t
175 size_i() const
176 {
177 return mSizeI ;
178 }
179
180 index_t
181 size_j() const
182 {
183 return mSizeJ ;
184 }
185
186 index_t
187 size_k() const
188 {
189 return mSizeK ;
190 }
191
192 index_t
193 size_l() const
194 {
195 return mSizeL ;
196 }
197 index_t
198 order() const
199 {
200 return mOrder ;
201 }
202
203//------------------------------------------------------------------------------
204// MEMORY
205//------------------------------------------------------------------------------
206
210 T *
211 data()
212 {
213 return mData ;
214 }
215
216// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
217
221 const T *
222 data() const
223 {
224 return mData ;
225 }
226
227//------------------------------------------------------------------------------
228// UTILITIES
229//------------------------------------------------------------------------------
230
234 void
235 fill( const T aValue )
236 {
237 // populate pointer array
238 std::fill_n( mData, mCapacity, aValue );
239 }
240
241//------------------------------------------------------------------------------
242
246 void
247 fill( const T aA, const T aB )
248 {
249 BELFEM_ASSERT( this->is_3333(), "Tensor must be 3x3x3x3" );
250 tensor::fill( mData, aA, aB );
251 }
252
253//------------------------------------------------------------------------------
254
258 void
259 fill_isotropic_elasticity( const T E, const T nu )
260 {
261 BELFEM_ASSERT( this->is_3333(), "Tensor must be 3x3x3x3" );
262
263 // bulk modulus
264 T K = E / ( 3. * ( 1.0 - 2.0 * nu ) );
265
266 // shear modulus
267 T G = E / ( 2.0 * ( 1.0 + nu ) );
268
269 tensor::fill( mData, K, G );
270 }
271
272//------------------------------------------------------------------------------
273
274 void
275 fill_orthotropic_elasticity(
276 const T aYoung1,
277 const T aYoung2,
278 const T aYoung3,
279 const T aPoisson23,
280 const T aPoisson13,
281 const T aPoisson12,
282 const T aShear23,
283 const T aShear13,
284 const T aShear12
285 )
286 {
287 BELFEM_ASSERT( this->is_3333(), "Tensor must be 3x3x3x3" );
288
289 // compliance matrix
290 Matrix< T > tS( 6, 6 );
291 compliance_matrix( aYoung1, aYoung2, aYoung3,
292 aPoisson23, aPoisson13, aPoisson12,
293 aShear23, aShear13, aShear12, tS );
294
295 // elasticity matrix
296 Matrix< T > tC = inv( tS );
297
298 // polulate data
299 tensor::mat_to_ten( tC, mData );
300 }
301
302//------------------------------------------------------------------------------
303
307 index_t
308 capacity() const
309 {
310 return mCapacity ;
311 }
312
313//------------------------------------------------------------------------------
314
318 bool
319 is_3333() const
320 {
321 return mSizeI == 3 && mSizeJ == 3 && mSizeK == 3 && mSizeL == 3 ;
322 }
323
324//------------------------------------------------------------------------------
325
326 void
327 print( const string aLabel="Tensor")
328 {
329 index_t tCount = 0 ;
330 fprintf( stdout, " %s:\n", aLabel.c_str() );
331
332 for( index_t l=0; l<mSizeL; ++l )
333 {
334 for( index_t k=0; k<mSizeK; ++k )
335 {
336 for( index_t j=0; j<mSizeJ; ++j )
337 {
338 for( index_t i=0; i<mSizeI; ++i )
339 {
340 fprintf( stdout, " %u : ( %u, %u, %u, %u ) = %12.3f\n",
341 ( unsigned int ) tCount,
342 ( unsigned int ) i,
343 ( unsigned int ) j,
344 ( unsigned int ) k,
345 ( unsigned int ) l,
346 ( double ) mData[ tCount ] );
347 ++tCount ;
348 }
349 }
350 }
351 }
352 }
353
354//------------------------------------------------------------------------------
355// ACCESS OPERATORS
356//------------------------------------------------------------------------------
357
361 T &
362 operator()(
363 const index_t I,
364 const index_t J,
365 const index_t K )
366 {
367 BELFEM_ASSERT( mOrder == 3, "Tensor order must be of order 3" );
368
369 BELFEM_ASSERT( I < mSizeI, "Index i out of bounds ( %u vs %u )",
370 ( unsigned int ) I,
371 ( unsigned int ) mSizeI );
372
373 BELFEM_ASSERT( J < mSizeJ, "Index j out of bounds ( %u vs %u )",
374 ( unsigned int ) J,
375 ( unsigned int ) mSizeJ );
376
377 BELFEM_ASSERT( K < mSizeK, "Index k out of bounds ( %u vs %u )",
378 ( unsigned int ) K,
379 ( unsigned int ) mSizeK );
380
381 return mData[ K * mOffsetK + J * mOffsetJ + I ] ;
382 }
383
384// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
385
389 const T &
390 operator()(
391 const index_t I,
392 const index_t J,
393 const index_t K ) const
394 {
395 BELFEM_ASSERT( mOrder == 3, "Tensor order must be of order 3" );
396
397 BELFEM_ASSERT( I < mSizeI, "Index i out of bounds ( %u vs %u )",
398 ( unsigned int ) I,
399 ( unsigned int ) mSizeI );
400
401 BELFEM_ASSERT( J < mSizeJ, "Index j out of bounds ( %u vs %u )",
402 ( unsigned int ) J,
403 ( unsigned int ) mSizeJ );
404
405 BELFEM_ASSERT( K < mSizeK, "Index k out of bounds ( %u vs %u )",
406 ( unsigned int ) K,
407 ( unsigned int ) mSizeK );
408
409 return mData[ K * mOffsetK + J * mOffsetJ + I ] ;
410 }
411
412// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
413
417 T &
418 operator()(
419 const index_t I,
420 const index_t J,
421 const index_t K,
422 const index_t L )
423 {
424 BELFEM_ASSERT( mOrder == 4, "Tensor order must be of order 4" );
425
426 BELFEM_ASSERT( I < mSizeI, "Index i out of bounds ( %u vs %u )",
427 ( unsigned int ) I,
428 ( unsigned int ) mSizeI );
429
430 BELFEM_ASSERT( J < mSizeJ, "Index j out of bounds ( %u vs %u )",
431 ( unsigned int ) J,
432 ( unsigned int ) mSizeJ );
433
434 BELFEM_ASSERT( K < mSizeK, "Index k out of bounds ( %u vs %u )",
435 ( unsigned int ) K,
436 ( unsigned int ) mSizeK );
437
438 BELFEM_ASSERT( L < mSizeL, "Index l out of bounds ( %u vs %u )",
439 ( unsigned int ) L,
440 ( unsigned int ) mSizeL );
441
442 return mData[ L * mOffsetL + K * mOffsetK + J * mOffsetJ + I ] ;
443 }
444
445// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
446
450 const T &
451 operator()(
452 const index_t I,
453 const index_t J,
454 const index_t K,
455 const index_t L ) const
456 {
457 BELFEM_ASSERT( mOrder == 4, "Tensor order must be of order 4" );
458
459 BELFEM_ASSERT( I < mSizeI, "Index i out of bounds ( %u vs %u )",
460 ( unsigned int ) I,
461 ( unsigned int ) mSizeI );
462
463 BELFEM_ASSERT( J < mSizeJ, "Index j out of bounds ( %u vs %u )",
464 ( unsigned int ) J,
465 ( unsigned int ) mSizeJ );
466
467 BELFEM_ASSERT( K < mSizeK, "Index k out of bounds ( %u vs %u )",
468 ( unsigned int ) K,
469 ( unsigned int ) mSizeK );
470
471 BELFEM_ASSERT( L < mSizeL, "Index l out of bounds ( %u vs %u )",
472 ( unsigned int ) L,
473 ( unsigned int ) mSizeL );
474
475 return mData[ L * mOffsetL + K * mOffsetK + J * mOffsetJ + I ] ;
476 }
477
478//------------------------------------------------------------------------------
479// Equal Operators
480//------------------------------------------------------------------------------
481
485 Tensor< T > &
486 operator=( const Tensor< T > & aTensor )
487 {
488 if( this == &aTensor ) return *this;
489
490 BELFEM_ASSERT( mOrder == aTensor.order(),
491 "Tensor orders must match for copy assignment" );
492
493 BELFEM_ASSERT( mSizeI == aTensor.size_i(),
494 "Tensor dimensions must match for copy assignment" );
495
496 BELFEM_ASSERT( mSizeJ == aTensor.size_j(),
497 "Tensor dimensions must match for copy assignment" );
498
499 BELFEM_ASSERT( mSizeK == aTensor.size_k(),
500 "Tensor dimensions must match for copy assignment" );
501
502 BELFEM_ASSERT( mSizeL == aTensor.size_l(),
503 "Tensor dimensions must match for copy assignment" );
504
505 std::memcpy( mData, aTensor.data(), mCapacity * sizeof( T ) );
506 return *this;
507 }
508
509// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
510
514 Tensor< T > &
515 operator=( Tensor< T > && aTensor )
516 {
517 if( this == &aTensor ) return *this;
518
519 BELFEM_ASSERT( mOrder == aTensor.order(),
520 "Tensor orders must match for move assignment" );
521
522 BELFEM_ASSERT( mSizeI == aTensor.size_i(),
523 "Tensor dimensions must match for move assignment" );
524
525 BELFEM_ASSERT( mSizeJ == aTensor.size_j(),
526 "Tensor dimensions must match for move assignment" );
527
528 BELFEM_ASSERT( mSizeK == aTensor.size_k(),
529 "Tensor dimensions must match for move assignment" );
530
531 BELFEM_ASSERT( mSizeL == aTensor.size_l(),
532 "Tensor dimensions must match for move assignment" );
533
534 free( mData );
535 mData = aTensor.mData ;
536 aTensor.mData = nullptr ;
537 return *this;
538 }
539
540// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
541
542 Tensor< T > &
543 operator=( const T & aScalar )
544 {
545 this->fill( aScalar );
546 return *this;
547 }
548
549// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
550
551 Tensor< T > &
552 operator=( const Matrix< T > & aElasticityMatrix )
553 {
554 BELFEM_ASSERT( this->is_3333(), "The tensor must be a 3x3x3x3 tensor" );
555
556 tensor::mat_to_ten( aElasticityMatrix, mData );
557
558 return *this;
559 }
560
561//------------------------------------------------------------------------------
562// Addition operators
563//------------------------------------------------------------------------------
564
565 Tensor< T > &
566 operator+=( const T & aScalar )
567 {
568 std::for_each( mData, mData + mCapacity,
569 [ aScalar ]( T & tVal )
570 { tVal += aScalar; } );
571 return *this;
572 }
573
574// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
575
576 Tensor< T > &
577 operator+=( const Tensor< T > & aTensor )
578 {
579 BELFEM_ASSERT( mOrder == aTensor.order(),
580 "Tensor orders must match for addition operator" );
581
582 BELFEM_ASSERT( mSizeI == aTensor.size_i(),
583 "Tensor dimensions must match for addition operator" );
584
585 BELFEM_ASSERT( mSizeJ == aTensor.size_j(),
586 "Tensor dimensions must match for addition operator" );
587
588 BELFEM_ASSERT( mSizeK == aTensor.size_k(),
589 "Tensor dimensions must match for addition operator" );
590
591 BELFEM_ASSERT( mSizeL == aTensor.size_l(),
592 "Tensor dimensions must match for addition operator" );
593
594 std::transform( mData, mData + mCapacity, aTensor.data(),
595 mData, std::plus< T >() );
596 return *this;
597 }
598
599//------------------------------------------------------------------------------
600// Subtraction operators
601//------------------------------------------------------------------------------
602
603 Tensor< T > &
604 operator-=( const T & aScalar )
605 {
606 std::for_each( mData, mData + mCapacity,
607 [ aScalar ]( T & tVal )
608 { tVal -= aScalar; } );
609 return *this;
610 }
611
612// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
613
614 Tensor< T > &
615 operator-=( const Tensor< T > & aTensor )
616 {
617 BELFEM_ASSERT( mOrder == aTensor.order(),
618 "Tensor orders must match for subtraction operator" );
619
620 BELFEM_ASSERT( mSizeI == aTensor.size_i(),
621 "Tensor dimensions must match for subtraction operator" );
622
623 BELFEM_ASSERT( mSizeJ == aTensor.size_j(),
624 "Tensor dimensions must match for subtraction operator" );
625
626 BELFEM_ASSERT( mSizeK == aTensor.size_k(),
627 "Tensor dimensions must match for subtraction operator" );
628
629 BELFEM_ASSERT( mSizeL == aTensor.size_l(),
630 "Tensor dimensions must match for subtraction operator" );
631
632 std::transform( mData, mData + mCapacity, aTensor.data(),
633 mData, std::minus< T >() );
634 return *this;
635 }
636
637//------------------------------------------------------------------------------
638// Multiplication
639//------------------------------------------------------------------------------
640
641 Tensor< T > &
642 operator*=( const T & aScalar )
643 {
644 std::for_each( mData, mData + mCapacity,
645 [ aScalar ]( T & tVal )
646 { tVal *= aScalar; } );
647 return *this;
648 }
649
650//------------------------------------------------------------------------------
651// Division
652//------------------------------------------------------------------------------
653
654 Tensor< T > &
655 operator/=( const T & aScalar )
656 {
657 std::for_each( mData, mData + mCapacity,
658 [ aScalar ]( T & tVal )
659 { tVal /= aScalar; } );
660 return *this;
661 }
662
663//------------------------------------------------------------------------------
664// Contraction
665//------------------------------------------------------------------------------
666
672 inline
673 void
674 ddot( const Tensor< T > & aB, Tensor< T > & aC )
675 {
676 BELFEM_ASSERT( this->is_3333(), "operating tensor must be 3x3x3x3" );
677 BELFEM_ASSERT( aB.is_3333(), "argument tensor must be 3x3x3x3" );
678 BELFEM_ASSERT( aC.is_3333(), "target tensor must be 3x3x3x3" );
679
680 tensor::contract44( mData, aB.data(), aC.data() );
681 }
682
683//------------------------------------------------------------------------------
684
690 inline
691 void
692 ddot( const Matrix< T > & aB, Matrix< T > & aC )
693 {
694 BELFEM_ASSERT( this->is_3333(),
695 "operating tensor must be 3x3x3x3" );
696
697 BELFEM_ASSERT( aB.n_rows() == 3
698 && aB.n_cols() == 3,
699 "argument matrix must be allocated as 3x3" );
700
701 BELFEM_ASSERT( aC.n_rows() == 3
702 && aC.n_cols() == 3,
703 "target matrix must be allocated as 3x3" );
704
705 tensor::contract42( mData, aB, aC );
706 }
707
708//------------------------------------------------------------------------------
709// Conversion
710//------------------------------------------------------------------------------
711
716 void
717 to_matrix( Matrix< real > & aMatrix )
718 {
719 BELFEM_ASSERT( aMatrix.n_rows() == 6
720 && aMatrix.n_cols() == 6,
721 "Matrix must be allocated as 6x6" );
722
723 BELFEM_ASSERT( this->is_3333(), "The tensor must be a 3x3x3x3 tensor" );
724
725 tensor::ten_to_mat( mData, aMatrix );
726 }
727 };
728//------------------------------------------------------------------------------
729
730 template< typename T >
731 Tensor< T >
732 operator+( const Tensor< T > & aB,
733 const Tensor< T > & aC )
734 {
735 BELFEM_ASSERT( aB.order() == aC.order(),
736 "Tensor orders must match for addition operator" );
737
738 BELFEM_ASSERT( aB.size_i() == aC.size_i(),
739 "Tensor dimensions must match for addition operator" );
740
741 BELFEM_ASSERT( aB.size_j() == aC.size_j(),
742 "Tensor dimensions must match for addition operator" );
743
744 BELFEM_ASSERT( aB.size_k() == aC.size_k(),
745 "Tensor dimensions must match for addition operator" );
746
747 BELFEM_ASSERT( aB.size_l() == aC.size_l(),
748 "Tensor dimensions must match for addition operator" );
749
750 Tensor< T > aA( aB );
751 aA += aC;
752
753 return aA ;
754 }
755
756//------------------------------------------------------------------------------
757
758 template< typename T >
759 Tensor< T >
760 operator-( const Tensor< T > & aB,
761 const Tensor< T > & aC )
762 {
763 BELFEM_ASSERT( aB.order() == aC.order(),
764 "Tensor orders must match for subtraction operator" );
765
766 BELFEM_ASSERT( aB.size_i() == aC.size_i(),
767 "Tensor dimensions must match for subtraction operator" );
768
769 BELFEM_ASSERT( aB.size_j() == aC.size_j(),
770 "Tensor dimensions must match for subtraction operator" );
771
772 BELFEM_ASSERT( aB.size_k() == aC.size_k(),
773 "Tensor dimensions must match for subtraction operator" );
774
775 BELFEM_ASSERT( aB.size_l() == aC.size_l(),
776 "Tensor dimensions must match for subtraction operator" );
777
778 Tensor< T > aA( aB );
779 aA -= aC;
780 return aA ;
781 }
782
783//------------------------------------------------------------------------------
784// contraction
785//------------------------------------------------------------------------------
786
787 template< typename T >
788 Tensor< T >
789 operator%( const Tensor< T > & aA,
790 const Tensor< T > & aB )
791 {
792 BELFEM_ASSERT( aA.is_3333() && aB.is_3333(),
793 "Both tensors must be of size 3x3x3x3" );
794
795 Tensor< T > aC( 3, 3, 3, 3 );
796 tensor::contract44( aA.data(), aB.data(), aC.data() );
797 return aC ;
798 }
799
800//------------------------------------------------------------------------------
801
802 template< typename T >
803 Matrix< T >
804 operator%( const Tensor< T > & aA,
805 const Matrix< T > & aB )
806 {
807 BELFEM_ASSERT( aA.is_3333(),
808 "Tensor A must be of size 3x3x3x3" );
809
810 BELFEM_ASSERT( aB.n_rows() == 3 && aB.n_cols() == 3,
811 "when contracting a 3x3x3x3 tensor with a matrix, latter one must be 3x3" );
812
813 Matrix< T > aC ( 3, 3 );
814 tensor::contract42( aA.data(), aB, aC );
815 return aC ;
816 }
817
818//------------------------------------------------------------------------------
819
820 template< typename T >
821 bool
822 operator==( const Tensor< T > & aA,
823 const Tensor< T > & aB )
824 {
825 if( aA.order() != aB.order()
826 || aA.size_i() != aB.size_i()
827 || aA.size_j() != aB.size_j()
828 || aA.size_k() != aB.size_k()
829 || aA.size_l() != aB.size_l() ) return false;
830
831 return tensor::equal_equal( aA.data(), aB.data(), aA.capacity() );
832 }
833
834//------------------------------------------------------------------------------
835} /* namespace belfem */
836#endif //BELFEM_CL_TENSOR_HPP
BELFEM_ASSERT
#define BELFEM_ASSERT(aCheck,...)
Definition assert.hpp:244
cl_Matrix.hpp
belfem::Matrix::n_rows
size_t n_rows() const
Definition cl_AR_Matrix.hpp:205
belfem::Matrix::n_cols
size_t n_cols() const
Definition cl_AR_Matrix.hpp:213
belfem::Tensor
Third- or fourth-order tensor container; the constitutive helpers (contraction, rotation,...
Definition cl_Tensor.hpp:43
belfem::Tensor::Tensor
Tensor(Tensor< T > &&aOther)
Definition cl_Tensor.hpp:148
belfem::Tensor::operator=
Tensor< T > & operator=(const T &aScalar)
Definition cl_Tensor.hpp:543
belfem::Tensor::Tensor
Tensor(const Tensor< T > &aOther)
Definition cl_Tensor.hpp:130
belfem::Tensor::operator+=
Tensor< T > & operator+=(const Tensor< T > &aTensor)
Definition cl_Tensor.hpp:577
belfem::Tensor::fill
void fill(const T aA, const T aB)
fill tensor in an isotropic way
Definition cl_Tensor.hpp:247
belfem::Tensor::operator+=
Tensor< T > & operator+=(const T &aScalar)
Definition cl_Tensor.hpp:566
belfem::Tensor::order
index_t order() const
Definition cl_Tensor.hpp:198
belfem::Tensor::to_matrix
void to_matrix(Matrix< real > &aMatrix)
converts a tensor to the elastitity matrix in Voigt notation
Definition cl_Tensor.hpp:717
belfem::Tensor::ddot
void ddot(const Matrix< T > &aB, Matrix< T > &aC)
contract with 3x3 matrox
Definition cl_Tensor.hpp:692
belfem::Tensor::operator*=
Tensor< T > & operator*=(const T &aScalar)
Definition cl_Tensor.hpp:642
belfem::Tensor::operator=
Tensor< T > & operator=(Tensor< T > &&aTensor)
move assignment
Definition cl_Tensor.hpp:515
belfem::Tensor::data
T * data()
expose the underlying raw pointer
Definition cl_Tensor.hpp:211
belfem::Tensor::operator=
Tensor< T > & operator=(const Matrix< T > &aElasticityMatrix)
Definition cl_Tensor.hpp:552
belfem::Tensor::size_j
index_t size_j() const
Definition cl_Tensor.hpp:181
belfem::Tensor::Tensor
Tensor(const index_t aSizeI, const index_t aSizeJ, const index_t aSizeK, const index_t aSizeL, const real aValue)
Definition cl_Tensor.hpp:106
belfem::Tensor::fill_orthotropic_elasticity
void fill_orthotropic_elasticity(const T aYoung1, const T aYoung2, const T aYoung3, const T aPoisson23, const T aPoisson13, const T aPoisson12, const T aShear23, const T aShear13, const T aShear12)
Definition cl_Tensor.hpp:275
belfem::Tensor::~Tensor
~Tensor()
Definition cl_Tensor.hpp:165
belfem::Tensor::ddot
void ddot(const Tensor< T > &aB, Tensor< T > &aC)
contract with other tensor
Definition cl_Tensor.hpp:674
belfem::Tensor::is_3333
bool is_3333() const
returns true if this is a 3x3x3x3 tensor
Definition cl_Tensor.hpp:319
belfem::Tensor::operator()
T & operator()(const index_t I, const index_t J, const index_t K)
access operator ( writable version )
Definition cl_Tensor.hpp:362
belfem::Tensor::print
void print(const string aLabel="Tensor")
Definition cl_Tensor.hpp:327
belfem::Tensor::size_l
index_t size_l() const
Definition cl_Tensor.hpp:193
belfem::Tensor::capacity
index_t capacity() const
memory size
Definition cl_Tensor.hpp:308
belfem::Tensor::operator=
Tensor< T > & operator=(const Tensor< T > &aTensor)
copy assignment
Definition cl_Tensor.hpp:486
belfem::Tensor::operator()
const T & operator()(const index_t I, const index_t J, const index_t K, const index_t L) const
access operator ( const version )
Definition cl_Tensor.hpp:451
belfem::Tensor::operator-=
Tensor< T > & operator-=(const T &aScalar)
Definition cl_Tensor.hpp:604
belfem::Tensor::operator-=
Tensor< T > & operator-=(const Tensor< T > &aTensor)
Definition cl_Tensor.hpp:615
belfem::Tensor::Tensor
Tensor(const index_t aSizeI, const index_t aSizeJ, const index_t aSizeK, const index_t aSizeL)
Definition cl_Tensor.hpp:67
belfem::Tensor::Tensor
Tensor(const Matrix< real > &aElasticityMatrix)
Definition cl_Tensor.hpp:117
belfem::Tensor::operator()
const T & operator()(const index_t I, const index_t J, const index_t K) const
access operator ( const version )
Definition cl_Tensor.hpp:390
belfem::Tensor::operator()
T & operator()(const index_t I, const index_t J, const index_t K, const index_t L)
access operator ( writable version )
Definition cl_Tensor.hpp:418
belfem::Tensor::fill
void fill(const T aValue)
fill all values
Definition cl_Tensor.hpp:235
belfem::Tensor::Tensor
Tensor(const index_t aSizeI, const index_t aSizeJ, const index_t aSizeK)
Definition cl_Tensor.hpp:87
belfem::Tensor::operator/=
Tensor< T > & operator/=(const T &aScalar)
Definition cl_Tensor.hpp:655
belfem::Tensor::fill_isotropic_elasticity
void fill_isotropic_elasticity(const T E, const T nu)
special funcition to create an isotropic elasticity tensor
Definition cl_Tensor.hpp:259
belfem::Tensor::size_k
index_t size_k() const
Definition cl_Tensor.hpp:187
belfem::Tensor::data
const T * data() const
expose the underlying raw pointer ( const version )
Definition cl_Tensor.hpp:222
belfem::Tensor::size_i
index_t size_i() const
Definition cl_Tensor.hpp:175
fn_TR_contract42.hpp
fn_TR_contract44.hpp
fn_TR_equal_equal.hpp
fn_TR_fill.hpp
fn_TR_mat_to_ten.hpp
fn_TR_ten_to_mat.hpp
fn_compliance_matrix.hpp
fn_inv.hpp
Inverse of a square matrix.
belfem::tensor::fill
void fill(T *A, const T a, const T b)
Definition fn_TR_fill.hpp:29
belfem::tensor::contract42
void contract42(const T *A, const T *B, T *C)
tensor contraction A_ijkl * B_kl = C_ij
Definition fn_TR_contract42_arma.hpp:28
belfem::tensor::mat_to_ten
void mat_to_ten(const T *C, T *A)
elasticity matrix conversion elasticity matrix C => tensor A
Definition fn_TR_mat_to_ten.hpp:27
belfem::tensor::contract44
void contract44(const T *A, const T *B, T *C)
tensor contraction A_ijmn * B_mnkl = C_ijkl
Definition fn_TR_contract44.hpp:26
belfem::tensor::equal_equal
bool equal_equal(const T *A, const T *B, const index_t aCapacity)
Definition fn_TR_equal_equal.hpp:23
belfem::tensor::ten_to_mat
void ten_to_mat(const T *A, T *C)
elasticity matrix conversion Elasticity tensor A => elasticity matrix C
Definition fn_TR_ten_to_mat.hpp:28
belfem
USER GUIDES:
Definition cl_Capacitor.cpp:16
belfem::compliance_matrix
void compliance_matrix(const T &aE1, const T &aE2, const T &aE3, const T &aNu23, const T &aNu13, const T &aNu12, const T &aG23, const T &aG31, const T &aG12, Matrix< T > &aS)
Definition fn_compliance_matrix.hpp:19
belfem::operator+
auto operator+(const Matrix< T > &aA, const Matrix< T > &aB) -> decltype(aA.matrix_data()+aB.matrix_data())
Definition op_MatrixPlus.hpp:23
belfem::MaterialProperty::nu
@ nu
Definition cl_Material.hpp:156
belfem::MaterialProperty::E
@ E
Definition cl_Material.hpp:155
belfem::operator-
auto operator-(const Matrix< T > &aA, const Matrix< T > &aB) -> decltype(aA.matrix_data() - aB.matrix_data())
Definition op_MatrixMinus.hpp:23
belfem::index_t
uint32_t index_t
Definition typedefs.hpp:52
belfem::inv
auto inv(const T &aExpression) -> decltype(arma::inv(aExpression))
Definition fn_AR_inv.hpp:22
belfem::operator%
auto operator%(const Vector< T > &aA, const Vector< T > &aB) -> decltype(aA.vector_data() % aB.vector_data())
Definition op_AR_VectorElementwiseMultiplication.hpp:22
belfem::real
double real
Definition typedefs.hpp:36
belfem::operator==
bool operator==(const Matrix< T > &aA, const Matrix< T > &aB)
Definition op_AR_MatrixEqualEqual.hpp:23