Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #pragma once
11 :
12 : #include "DenseMatrix.h"
13 : #include "MooseError.h"
14 : #include "ADReal.h"
15 : #include "MooseTypes.h"
16 :
17 : #include "libmesh/type_tensor.h"
18 : #include "libmesh/dense_vector.h"
19 :
20 : #include "metaphysicl/raw_type.h"
21 :
22 : // C++ includes
23 : #include <iomanip>
24 :
25 : /**
26 : * This class defines a Tensor that can change its shape. This means
27 : * a 3x3x3x3 Tensor can be represented as a 9x9 or an 81x1. Further,
28 : * the values of this tensor are _COLUMN_ major ordered!
29 : */
30 : template <typename T>
31 : class ColumnMajorMatrixTempl
32 : {
33 : public:
34 : /**
35 : * Constructor that sets an initial number of entries and shape.
36 : * Defaults to creating the same size tensor as TensorValue
37 : */
38 : explicit ColumnMajorMatrixTempl(const unsigned int rows = Moose::dim,
39 : const unsigned int cols = Moose::dim);
40 :
41 : /**
42 : * Copy Constructor defined in terms of operator=()
43 : */
44 : ColumnMajorMatrixTempl(const ColumnMajorMatrixTempl<T> & rhs);
45 :
46 : /**
47 : * Constructor that fills in the ColumnMajorMatrixTempl with values from a libMesh TypeTensor
48 : */
49 : explicit ColumnMajorMatrixTempl(const TypeTensor<T> & tensor);
50 :
51 : explicit ColumnMajorMatrixTempl(const DenseMatrix<T> & rhs);
52 :
53 : explicit ColumnMajorMatrixTempl(const DenseVector<T> & rhs);
54 :
55 : /**
56 : * Constructor that takes in 3 vectors and uses them to create columns
57 : */
58 : ColumnMajorMatrixTempl(const TypeVector<T> & col1,
59 : const TypeVector<T> & col2,
60 : const TypeVector<T> & col3);
61 :
62 : /**
63 : * The total number of entries in the Tensor.
64 : * i.e. cols * rows
65 : */
66 : unsigned int numEntries() const;
67 :
68 : /**
69 : * Change the shape of the tensor.
70 : * Note that cols * rows should be equal to numEntries()!
71 : */
72 : void reshape(const unsigned int rows, const unsigned int cols);
73 :
74 : /**
75 : * Get the i,j entry
76 : * j defaults to zero so you can use it as a column vector.
77 : */
78 : T & operator()(const unsigned int i, const unsigned int j = 0);
79 :
80 : /**
81 : * Get the i,j entry
82 : *
83 : * j defaults to zero so you can use it as a column vector.
84 : * This is the version used for a const ColumnMajorMatrixTempl.
85 : */
86 : T operator()(const unsigned int i, const unsigned int j = 0) const;
87 :
88 : /**
89 : * Print the tensor
90 : */
91 : void print();
92 :
93 : /**
94 : * Prints to file
95 : */
96 : void print_scientific(std::ostream & os);
97 :
98 : /**
99 : * Fills the passed in tensor with the values from this tensor.
100 : */
101 : void fill(TypeTensor<T> & tensor);
102 :
103 : /**
104 : * Fills the passed in dense matrix with the values from this tensor.
105 : */
106 : void fill(DenseMatrix<T> & rhs);
107 :
108 : /**
109 : * Fills the passed in dense vector with the values from this tensor.
110 : */
111 : void fill(DenseVector<T> & rhs);
112 :
113 : /**
114 : * Returns a matrix that is the transpose of the matrix this
115 : * was called on.
116 : */
117 : ColumnMajorMatrixTempl<T> transpose() const;
118 :
119 : /**
120 : * Returns a matrix that is the deviatoric of the matrix this
121 : * was called on.
122 : */
123 : ColumnMajorMatrixTempl<T> deviatoric();
124 :
125 : /**
126 : * Returns a matrix that is the absolute value of the matrix this
127 : * was called on.
128 : */
129 : ColumnMajorMatrixTempl<T> abs();
130 :
131 : /**
132 : * Set the value of each of the diagonals to the passed in value.
133 : */
134 : void setDiag(T value);
135 :
136 : /**
137 : * Add to each of the diagonals the passsed in value.
138 : */
139 : void addDiag(T value);
140 :
141 : /**
142 : * The trace of the CMM.
143 : */
144 : T tr() const;
145 :
146 : /**
147 : * Zero the matrix.
148 : */
149 : void zero();
150 :
151 : /**
152 : * Turn the matrix into an identity matrix.
153 : */
154 : void identity();
155 :
156 : /**
157 : * Double contraction of two matrices ie A : B = Sum(A_ab * B_ba)
158 : */
159 : T doubleContraction(const ColumnMajorMatrixTempl<T> & rhs) const;
160 :
161 : /**
162 : * The Euclidean norm of the matrix.
163 : */
164 : T norm();
165 :
166 : /**
167 : * Returns the number of rows
168 : */
169 : unsigned int n() const;
170 :
171 : /**
172 : * Returns the number of columns
173 : */
174 : unsigned int m() const;
175 :
176 : /**
177 : * Returns eigen system solve for a symmetric real matrix
178 : */
179 : void eigen(ColumnMajorMatrixTempl<T> & eval, ColumnMajorMatrixTempl<T> & evec) const;
180 :
181 : /**
182 : * Returns eigen system solve for a non-symmetric real matrix
183 : */
184 : void eigenNonsym(ColumnMajorMatrixTempl<T> & eval_real,
185 : ColumnMajorMatrixTempl<T> & eval_img,
186 : ColumnMajorMatrixTempl<T> & evec_right,
187 : ColumnMajorMatrixTempl<T> & eve_left) const;
188 :
189 : /**
190 : * Returns matrix that is the exponential of the matrix this was called on
191 : */
192 : void exp(ColumnMajorMatrixTempl<T> & z) const;
193 :
194 : /**
195 : * Returns inverse of a general matrix
196 : */
197 : void inverse(ColumnMajorMatrixTempl<T> & invA) const;
198 :
199 : /**
200 : * Returns a reference to the raw data pointer
201 : */
202 : T * rawData();
203 : const T * rawData() const;
204 :
205 : /**
206 : * Kronecker Product
207 : */
208 : ColumnMajorMatrixTempl<T> kronecker(const ColumnMajorMatrixTempl<T> & rhs) const;
209 :
210 : /**
211 : * Sets the values in _this_ tensor to the values on the RHS.
212 : * Will also reshape this tensor if necessary.
213 : */
214 : ColumnMajorMatrixTempl<T> & operator=(const TypeTensor<T> & rhs);
215 :
216 : /**
217 : * Sets the values in _this_ dense matrix to the values on the RHS.
218 : * Will also reshape this tensor if necessary.
219 : */
220 : ColumnMajorMatrixTempl<T> & operator=(const DenseMatrix<T> & rhs);
221 :
222 : /**
223 : * Sets the values in _this_ dense vector to the values on the RHS.
224 : * Will also reshape this tensor if necessary.
225 : */
226 : ColumnMajorMatrixTempl<T> & operator=(const DenseVector<T> & rhs);
227 :
228 : /**
229 : * Sets the values in _this_ tensor to the values on the RHS
230 : * Will also reshape this tensor if necessary.
231 : */
232 : template <typename T2>
233 : ColumnMajorMatrixTempl<T> & operator=(const ColumnMajorMatrixTempl<T2> & rhs);
234 :
235 : /**
236 : * defaulted operator=
237 : */
238 182 : ColumnMajorMatrixTempl<T> & operator=(const ColumnMajorMatrixTempl<T> & rhs) = default;
239 :
240 : /**
241 : * Scalar multiplication of the ColumnMajorMatrixTempl
242 : */
243 : ColumnMajorMatrixTempl<T> operator*(T scalar) const;
244 :
245 : /**
246 : * Matrix Vector Multiplication of the libMesh TypeVector Type
247 : */
248 : ColumnMajorMatrixTempl<T> operator*(const TypeVector<T> & rhs) const;
249 :
250 : // /**
251 : // * Matrix Vector Multiplication of the TypeTensor Product. Note that the
252 : // * Tensor type is treated as a single dimension Vector for this operation
253 : // */
254 : // ColumnMajorMatrixTempl operator*(const TypeTensor<T> & rhs) const;
255 :
256 : /**
257 : * Matrix Matrix Multiplication
258 : */
259 : ColumnMajorMatrixTempl<T> operator*(const ColumnMajorMatrixTempl<T> & rhs) const;
260 :
261 : /**
262 : * Matrix Matrix Addition
263 : */
264 : ColumnMajorMatrixTempl<T> operator+(const ColumnMajorMatrixTempl<T> & rhs) const;
265 :
266 : /**
267 : * Matrix Matrix Subtraction
268 : */
269 : ColumnMajorMatrixTempl<T> operator-(const ColumnMajorMatrixTempl<T> & rhs) const;
270 :
271 : /**
272 : * Matrix Matrix Addition plus assignment
273 : *
274 : * Note that this is faster than regular addition
275 : * because the result doesn't have to get copied out
276 : */
277 : ColumnMajorMatrixTempl<T> & operator+=(const ColumnMajorMatrixTempl<T> & rhs);
278 :
279 : /**
280 : * Matrix Tensor Addition Plus Assignment
281 : */
282 : ColumnMajorMatrixTempl<T> & operator+=(const TypeTensor<T> & rhs);
283 :
284 : /**
285 : * Matrix Matrix Subtraction plus assignment
286 : *
287 : * Note that this is faster than regular subtraction
288 : * because the result doesn't have to get copied out
289 : */
290 : ColumnMajorMatrixTempl<T> & operator-=(const ColumnMajorMatrixTempl<T> & rhs);
291 :
292 : /**
293 : * Scalar addition
294 : */
295 : ColumnMajorMatrixTempl<T> operator+(T scalar) const;
296 :
297 : /**
298 : * Scalar Multiplication plus assignment
299 : */
300 : ColumnMajorMatrixTempl<T> & operator*=(T scalar);
301 :
302 : /**
303 : * Scalar Division plus assignment
304 : */
305 : ColumnMajorMatrixTempl<T> & operator/=(T scalar);
306 :
307 : /**
308 : * Scalar Addition plus assignment
309 : */
310 : ColumnMajorMatrixTempl<T> & operator+=(T scalar);
311 :
312 : /**
313 : * Check if matrix is square
314 : */
315 : void checkSquareness() const;
316 :
317 : /**
318 : * Check if matrices are of same shape
319 : */
320 : void checkShapeEquality(const ColumnMajorMatrixTempl<T> & rhs) const;
321 :
322 : /**
323 : * Equality operators
324 : */
325 : bool operator==(const ColumnMajorMatrixTempl<T> & rhs) const;
326 : bool operator!=(const ColumnMajorMatrixTempl<T> & rhs) const;
327 :
328 : protected:
329 : unsigned int _n_rows, _n_cols, _n_entries;
330 : std::vector<T> _values;
331 :
332 : template <typename T2>
333 : friend void dataStore(std::ostream &, ColumnMajorMatrixTempl<T2> &, void *);
334 : template <typename T2>
335 : friend void dataLoad(std::istream &, ColumnMajorMatrixTempl<T2> &, void *);
336 : };
337 :
338 : template <typename T>
339 : inline unsigned int
340 7 : ColumnMajorMatrixTempl<T>::numEntries() const
341 : {
342 7 : return _n_entries;
343 : }
344 :
345 : template <typename T>
346 : inline void
347 159 : ColumnMajorMatrixTempl<T>::reshape(unsigned int rows, unsigned int cols)
348 : {
349 159 : if (cols * rows == _n_entries)
350 : {
351 159 : _n_rows = rows;
352 159 : _n_cols = cols;
353 : }
354 : else
355 : {
356 0 : _n_rows = rows;
357 0 : _n_cols = cols;
358 0 : _n_entries = _n_rows * _n_cols;
359 0 : _values.resize(_n_entries);
360 : }
361 159 : }
362 :
363 : template <typename T>
364 : inline T &
365 100013957 : ColumnMajorMatrixTempl<T>::operator()(const unsigned int i, const unsigned int j)
366 : {
367 100013957 : if ((i * j) >= _n_entries)
368 1 : mooseError("Reference outside of ColumnMajorMatrix bounds!");
369 :
370 : // Row major indexing!
371 100013956 : return _values[(j * _n_rows) + i];
372 : }
373 :
374 : template <typename T>
375 : inline T
376 12934 : ColumnMajorMatrixTempl<T>::operator()(const unsigned int i, const unsigned int j) const
377 : {
378 12934 : if ((i * j) >= _n_entries)
379 0 : mooseError("Reference outside of ColumnMajorMatrix bounds!");
380 :
381 : // Row major indexing!
382 12934 : return _values[(j * _n_rows) + i];
383 : }
384 :
385 : template <typename T>
386 : inline void
387 0 : ColumnMajorMatrixTempl<T>::print()
388 : {
389 0 : ColumnMajorMatrixTempl<T> & s = (*this);
390 :
391 0 : for (unsigned int i = 0; i < _n_rows; i++)
392 : {
393 0 : for (unsigned int j = 0; j < _n_cols; j++)
394 0 : Moose::out << std::setw(15) << s(i, j) << " ";
395 :
396 0 : Moose::out << std::endl;
397 : }
398 0 : }
399 :
400 : template <typename T>
401 : inline void
402 0 : ColumnMajorMatrixTempl<T>::print_scientific(std::ostream & os)
403 : {
404 0 : ColumnMajorMatrixTempl<T> & s = (*this);
405 :
406 0 : for (unsigned int i = 0; i < _n_rows; i++)
407 : {
408 0 : for (unsigned int j = 0; j < _n_cols; j++)
409 0 : os << std::setw(15) << std::scientific << std::setprecision(8) << s(i, j) << " ";
410 :
411 0 : os << std::endl;
412 : }
413 0 : }
414 :
415 : template <typename T>
416 : inline void
417 2 : ColumnMajorMatrixTempl<T>::fill(TypeTensor<T> & tensor)
418 : {
419 2 : if (Moose::dim * Moose::dim != _n_entries)
420 1 : mooseError(
421 : "Cannot fill tensor! The ColumnMajorMatrix doesn't have the same number of entries!");
422 :
423 4 : for (const auto j : libMesh::make_range(Moose::dim))
424 12 : for (const auto i : libMesh::make_range(Moose::dim))
425 9 : tensor(i, j) = _values[j * Moose::dim + i];
426 1 : }
427 :
428 : template <typename T>
429 : inline void
430 1 : ColumnMajorMatrixTempl<T>::fill(DenseMatrix<T> & rhs)
431 : {
432 1 : if (rhs.n() * rhs.m() != _n_entries)
433 1 : mooseError(
434 : "Cannot fill dense matrix! The ColumnMajorMatrix doesn't have the same number of entries!");
435 :
436 0 : for (unsigned int j = 0, index = 0; j < rhs.m(); ++j)
437 0 : for (unsigned int i = 0; i < rhs.n(); ++i, ++index)
438 0 : rhs(i, j) = _values[index];
439 0 : }
440 :
441 : template <typename T>
442 : inline void
443 1 : ColumnMajorMatrixTempl<T>::fill(DenseVector<T> & rhs)
444 : {
445 1 : if (_n_rows != rhs.size() || _n_cols != 1)
446 1 : mooseError("ColumnMajorMatrix and DenseVector must be the same shape for a fill!");
447 :
448 0 : for (unsigned int i = 0; i < _n_rows; ++i)
449 0 : rhs(i) = (*this)(i);
450 0 : }
451 :
452 : template <typename T>
453 : inline ColumnMajorMatrixTempl<T>
454 1 : ColumnMajorMatrixTempl<T>::transpose() const
455 : {
456 1 : const ColumnMajorMatrixTempl<T> & s = (*this);
457 :
458 1 : ColumnMajorMatrixTempl<T> ret_matrix(_n_cols, _n_rows);
459 :
460 3 : for (unsigned int i = 0; i < _n_rows; i++)
461 6 : for (unsigned int j = 0; j < _n_cols; j++)
462 4 : ret_matrix(j, i) = s(i, j);
463 :
464 1 : return ret_matrix;
465 0 : }
466 :
467 : template <typename T>
468 : inline ColumnMajorMatrixTempl<T>
469 0 : ColumnMajorMatrixTempl<T>::deviatoric()
470 : {
471 0 : ColumnMajorMatrixTempl<T> & s = (*this);
472 :
473 0 : ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols), I(_n_rows, _n_cols);
474 :
475 0 : I.identity();
476 :
477 0 : for (unsigned int i = 0; i < _n_rows; i++)
478 0 : for (unsigned int j = 0; j < _n_cols; j++)
479 0 : ret_matrix(i, j) = s(i, j) - I(i, j) * (s.tr() / 3.0);
480 :
481 0 : return ret_matrix;
482 0 : }
483 :
484 : template <typename T>
485 : inline void
486 2 : ColumnMajorMatrixTempl<T>::setDiag(T value)
487 : {
488 2 : this->checkSquareness();
489 :
490 3 : for (unsigned int i = 0; i < _n_rows; i++)
491 2 : (*this)(i, i) = value;
492 1 : }
493 :
494 : template <typename T>
495 : inline void
496 0 : ColumnMajorMatrixTempl<T>::addDiag(T value)
497 : {
498 0 : this->checkSquareness();
499 :
500 0 : for (unsigned int i = 0; i < _n_rows; i++)
501 0 : (*this)(i, i) += value;
502 0 : }
503 :
504 : template <typename T>
505 : inline T
506 1 : ColumnMajorMatrixTempl<T>::tr() const
507 : {
508 1 : this->checkSquareness();
509 :
510 1 : T trace = 0;
511 :
512 3 : for (unsigned int i = 0; i < _n_rows; i++)
513 2 : trace += (*this)(i, i);
514 :
515 1 : return trace;
516 0 : }
517 :
518 : template <typename T>
519 : inline void
520 2 : ColumnMajorMatrixTempl<T>::zero()
521 : {
522 10 : for (unsigned int i = 0; i < _n_entries; i++)
523 8 : _values[i] = 0;
524 2 : }
525 :
526 : template <typename T>
527 : inline void
528 1 : ColumnMajorMatrixTempl<T>::identity()
529 : {
530 1 : this->checkSquareness();
531 :
532 1 : zero();
533 :
534 3 : for (unsigned int i = 0; i < _n_rows; i++)
535 2 : (*this)(i, i) = 1;
536 1 : }
537 :
538 : template <typename T>
539 : inline T
540 3 : ColumnMajorMatrixTempl<T>::doubleContraction(const ColumnMajorMatrixTempl<T> & rhs) const
541 : {
542 3 : this->checkShapeEquality(rhs);
543 :
544 2 : T value = 0;
545 :
546 4 : for (unsigned int j = 0; j < _n_cols; j++)
547 10 : for (unsigned int i = 0; i < _n_rows; i++)
548 8 : value += (*this)(i, j) * rhs(i, j);
549 :
550 2 : return value;
551 0 : }
552 :
553 : template <typename T>
554 : inline unsigned int
555 0 : ColumnMajorMatrixTempl<T>::n() const
556 : {
557 0 : return _n_rows;
558 : }
559 :
560 : template <typename T>
561 : inline unsigned int
562 0 : ColumnMajorMatrixTempl<T>::m() const
563 : {
564 0 : return _n_cols;
565 : }
566 :
567 : template <typename T>
568 : inline T *
569 20 : ColumnMajorMatrixTempl<T>::rawData()
570 : {
571 20 : return &_values[0];
572 : }
573 :
574 : template <typename T>
575 : inline const T *
576 0 : ColumnMajorMatrixTempl<T>::rawData() const
577 : {
578 0 : return &_values[0];
579 : }
580 :
581 : template <typename T>
582 : inline ColumnMajorMatrixTempl<T> &
583 1 : ColumnMajorMatrixTempl<T>::operator=(const TypeTensor<T> & rhs)
584 : {
585 : // Resize the tensor if necessary
586 1 : if ((Moose::dim * Moose::dim) != _n_entries)
587 : {
588 0 : _n_entries = Moose::dim * Moose::dim;
589 0 : _values.resize(_n_entries);
590 : }
591 :
592 : // Make sure the shape is correct
593 1 : reshape(Moose::dim, Moose::dim);
594 :
595 1 : ColumnMajorMatrixTempl<T> & s = (*this);
596 :
597 : // Copy the values
598 4 : for (unsigned int j = 0; j < _n_cols; j++)
599 12 : for (unsigned int i = 0; i < _n_cols; i++)
600 9 : s(i, j) = rhs(i, j);
601 :
602 1 : return *this;
603 : }
604 :
605 : template <typename T>
606 : template <typename T2>
607 : inline ColumnMajorMatrixTempl<T> &
608 : ColumnMajorMatrixTempl<T>::operator=(const ColumnMajorMatrixTempl<T2> & rhs)
609 : {
610 : this->reshape(rhs.m(), rhs.n());
611 :
612 : for (MooseIndex(rhs.m()) i = 0; i < rhs.m(); ++i)
613 : for (MooseIndex(rhs.n()) j = 0; j < rhs.n(); ++j)
614 : (*this)(i, j) = rhs(i, j);
615 :
616 : return *this;
617 : }
618 :
619 : template <typename T>
620 : inline ColumnMajorMatrixTempl<T>
621 1 : ColumnMajorMatrixTempl<T>::operator*(T scalar) const
622 : {
623 1 : ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols);
624 :
625 10 : for (unsigned int i = 0; i < _n_entries; i++)
626 9 : ret_matrix._values[i] = _values[i] * scalar;
627 :
628 1 : return ret_matrix;
629 0 : }
630 :
631 : template <typename T>
632 : inline ColumnMajorMatrixTempl<T>
633 1 : ColumnMajorMatrixTempl<T>::operator*(const TypeVector<T> & rhs) const
634 : {
635 1 : if (_n_cols != Moose::dim)
636 1 : mooseError("Cannot perform matvec operation! The column dimension of "
637 : "the ColumnMajorMatrix does not match the TypeVector!");
638 :
639 0 : ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, 1);
640 :
641 0 : for (unsigned int i = 0; i < _n_rows; ++i)
642 0 : for (unsigned int j = 0; j < _n_cols; ++j)
643 0 : ret_matrix._values[i] += (*this)(i, j) * rhs(j);
644 :
645 0 : return ret_matrix;
646 0 : }
647 :
648 : // template <typename T>
649 : // inline ColumnMajorMatrixTempl<T>
650 : // ColumnMajorMatrixTempl<T>::operator*(const TypeTensor<T> & rhs) const
651 : // {
652 : // mooseAssert(_n_cols == LIBMESH_DIM*LIBMESH_DIM, "Cannot perform matvec operation! The
653 : // ColumnMajorMatrixTempl<T> doesn't have the correct shape!");
654 :
655 : // ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, 1);
656 :
657 : // for (unsigned int i=0; i<_n_rows; ++i)
658 : // for (unsigned int j=0; j<_n_cols; ++j)
659 : // // Treat the Tensor as a column major column vector
660 : // ret_matrix._values[i] += (*this)(i, j) * rhs(j%3, j/3);
661 :
662 : // return ret_matrix;
663 : // }
664 :
665 : template <typename T>
666 : inline ColumnMajorMatrixTempl<T>
667 7 : ColumnMajorMatrixTempl<T>::operator*(const ColumnMajorMatrixTempl<T> & rhs) const
668 : {
669 7 : if (_n_cols != rhs._n_rows)
670 1 : mooseError(
671 : "Cannot perform matrix multiply! The shapes of the two operands are not compatible!");
672 :
673 6 : ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, rhs._n_cols);
674 :
675 24 : for (unsigned int i = 0; i < ret_matrix._n_rows; ++i)
676 66 : for (unsigned int j = 0; j < ret_matrix._n_cols; ++j)
677 192 : for (unsigned int k = 0; k < _n_cols; ++k)
678 144 : ret_matrix(i, j) += (*this)(i, k) * rhs(k, j);
679 :
680 6 : return ret_matrix;
681 0 : }
682 :
683 : template <typename T>
684 : inline ColumnMajorMatrixTempl<T>
685 1 : ColumnMajorMatrixTempl<T>::operator+(const ColumnMajorMatrixTempl<T> & rhs) const
686 : {
687 1 : this->checkShapeEquality(rhs);
688 :
689 1 : ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols);
690 :
691 7 : for (unsigned int i = 0; i < _n_entries; i++)
692 6 : ret_matrix._values[i] = _values[i] + rhs._values[i];
693 :
694 1 : return ret_matrix;
695 0 : }
696 :
697 : template <typename T>
698 : inline ColumnMajorMatrixTempl<T>
699 0 : ColumnMajorMatrixTempl<T>::operator-(const ColumnMajorMatrixTempl<T> & rhs) const
700 : {
701 0 : this->checkShapeEquality(rhs);
702 :
703 0 : ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols);
704 :
705 0 : for (unsigned int i = 0; i < _n_entries; i++)
706 0 : ret_matrix._values[i] = _values[i] - rhs._values[i];
707 :
708 0 : return ret_matrix;
709 0 : }
710 :
711 : template <typename T>
712 : inline ColumnMajorMatrixTempl<T> &
713 1 : ColumnMajorMatrixTempl<T>::operator+=(const ColumnMajorMatrixTempl<T> & rhs)
714 : {
715 1 : this->checkShapeEquality(rhs);
716 :
717 7 : for (unsigned int i = 0; i < _n_entries; i++)
718 6 : _values[i] += rhs._values[i];
719 :
720 1 : return *this;
721 : }
722 :
723 : template <typename T>
724 : inline ColumnMajorMatrixTempl<T> &
725 1 : ColumnMajorMatrixTempl<T>::operator+=(const TypeTensor<T> & rhs)
726 : {
727 1 : if ((_n_rows != Moose::dim) || (_n_cols != Moose::dim))
728 1 : mooseError("Cannot perform matrix addition and assignment! The shapes of the two operands are "
729 : "not compatible!");
730 :
731 0 : for (const auto j : libMesh::make_range(Moose::dim))
732 0 : for (const auto i : libMesh::make_range(Moose::dim))
733 0 : (*this)(i, j) += rhs(i, j);
734 :
735 0 : return *this;
736 : }
737 :
738 : template <typename T>
739 : inline ColumnMajorMatrixTempl<T> &
740 1 : ColumnMajorMatrixTempl<T>::operator-=(const ColumnMajorMatrixTempl<T> & rhs)
741 : {
742 1 : this->checkShapeEquality(rhs);
743 :
744 7 : for (unsigned int i = 0; i < _n_entries; i++)
745 6 : _values[i] -= rhs._values[i];
746 :
747 1 : return *this;
748 : }
749 :
750 : template <typename T>
751 : inline ColumnMajorMatrixTempl<T>
752 1 : ColumnMajorMatrixTempl<T>::operator+(T scalar) const
753 : {
754 1 : ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols);
755 :
756 10 : for (unsigned int i = 0; i < _n_entries; i++)
757 9 : ret_matrix._values[i] = _values[i] + scalar;
758 :
759 1 : return ret_matrix;
760 0 : }
761 :
762 : template <typename T>
763 : inline ColumnMajorMatrixTempl<T> &
764 1 : ColumnMajorMatrixTempl<T>::operator*=(T scalar)
765 : {
766 10 : for (unsigned int i = 0; i < _n_entries; i++)
767 9 : _values[i] *= scalar;
768 1 : return *this;
769 : }
770 :
771 : template <typename T>
772 : inline ColumnMajorMatrixTempl<T> &
773 1 : ColumnMajorMatrixTempl<T>::operator/=(T scalar)
774 : {
775 5 : for (unsigned int i = 0; i < _n_entries; i++)
776 4 : _values[i] /= scalar;
777 1 : return *this;
778 : }
779 :
780 : template <typename T>
781 : inline ColumnMajorMatrixTempl<T> &
782 1 : ColumnMajorMatrixTempl<T>::operator+=(T scalar)
783 : {
784 10 : for (unsigned int i = 0; i < _n_entries; i++)
785 9 : _values[i] += scalar;
786 1 : return *this;
787 : }
788 :
789 : template <typename T>
790 : inline bool
791 16 : ColumnMajorMatrixTempl<T>::operator==(const ColumnMajorMatrixTempl<T> & rhs) const
792 : {
793 16 : if (_n_entries != rhs._n_entries || _n_rows != rhs._n_rows || _n_cols != rhs._n_cols)
794 2 : return false;
795 14 : return std::equal(_values.begin(), _values.end(), rhs._values.begin());
796 : }
797 :
798 : template <typename T>
799 : inline bool
800 2 : ColumnMajorMatrixTempl<T>::operator!=(const ColumnMajorMatrixTempl<T> & rhs) const
801 : {
802 2 : return !(*this == rhs);
803 : }
804 :
805 : typedef ColumnMajorMatrixTempl<Real> ColumnMajorMatrix;
806 :
807 : namespace MetaPhysicL
808 : {
809 : template <typename T>
810 : struct RawType<ColumnMajorMatrixTempl<T>>
811 : {
812 : typedef ColumnMajorMatrixTempl<typename RawType<T>::value_type> value_type;
813 :
814 : static value_type value(const ColumnMajorMatrixTempl<T> & in)
815 : {
816 : value_type ret(in.m(), in.n());
817 : for (MooseIndex(in.m()) i = 0; i < in.m(); ++i)
818 : for (MooseIndex(in.n()) j = 0; j < in.n(); ++j)
819 : ret(i, j) = raw_value(in(i, j));
820 :
821 : return ret;
822 : }
823 : };
824 : }
|
