Line data Source code
1 : #ifndef SET_MATRIX4_DIAGONAL_H
2 : #define SET_MATRIX4_DIAGONAL_H
3 :
4 : #include "AMReX_Config.H"
5 : #include "Util/Util.H"
6 : #include "Base.H"
7 :
8 : #include "Set.H"
9 : #include "Set/Matrix3.H"
10 : #include "Set/Matrix4.H"
11 :
12 : namespace Set
13 : {
14 : template<>
15 : class Matrix4<AMREX_SPACEDIM, Sym::Diagonal>
16 : {
17 : std::array<std::array<Set::Scalar,AMREX_SPACEDIM>,AMREX_SPACEDIM> A = {{
18 : AMREX_D_DECL(
19 : {AMREX_D_DECL(NAN,NAN,NAN)},
20 : {AMREX_D_DECL(NAN,NAN,NAN)},
21 : {AMREX_D_DECL(NAN,NAN,NAN)}
22 : )}};
23 :
24 :
25 : public:
26 0 : Set::Matrix getA() const
27 : {
28 0 : Set::Matrix ret;
29 0 : for (int i = 0 ; i < AMREX_SPACEDIM; i++)
30 0 : for (int j = 0 ; j < AMREX_SPACEDIM; j++)
31 0 : ret(i,j)=A[i][j];
32 0 : return ret;
33 : }
34 :
35 110811 : AMREX_GPU_HOST_DEVICE Matrix4() = default;
36 : //AMREX_GPU_HOST_DEVICE Matrix4(Set::Matrix a_A): A(a_A) {};
37 : AMREX_FORCE_INLINE
38 : Scalar operator () (const int i, const int j, const int k, const int l) const
39 : {
40 1940 : if (i == k && j == l) return A[i][j];
41 1680 : else return 0.0;
42 : }
43 44 : void Randomize()
44 : {
45 154 : for (int i = 0 ; i < AMREX_SPACEDIM; i++)
46 396 : for (int j = 0 ; j < AMREX_SPACEDIM; j++)
47 286 : A[i][j] = Util::Random();
48 44 : }
49 2 : static Matrix4<AMREX_SPACEDIM,Sym::Diagonal> Random()
50 : {
51 2 : Matrix4<AMREX_SPACEDIM,Sym::Diagonal> ret;
52 2 : ret.Randomize();
53 2 : return ret;
54 : }
55 0 : void Print(std::ostream& os)
56 : {
57 0 : os << "A = \n" << getA();
58 0 : }
59 27649 : static Matrix4<AMREX_SPACEDIM, Sym::Diagonal> Zero()
60 : {
61 27649 : Matrix4<AMREX_SPACEDIM, Sym::Diagonal> zero;
62 82958 : for (int i = 0 ; i < AMREX_SPACEDIM; i++)
63 165960 : for (int j = 0 ; j < AMREX_SPACEDIM; j++)
64 110651 : zero.A[i][j] = 0.0;
65 27649 : return zero;
66 : }
67 9 : static Matrix4<AMREX_SPACEDIM, Sym::Diagonal> Ones()
68 : {
69 9 : Matrix4<AMREX_SPACEDIM, Sym::Diagonal> ones;
70 27 : for (int i = 0 ; i < AMREX_SPACEDIM; i++)
71 54 : for (int j = 0 ; j < AMREX_SPACEDIM; j++)
72 36 : ones.A[i][j] = 1.0;
73 9 : return ones;
74 : }
75 : friend Set::Matrix operator * (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a, const Set::Matrix& b);
76 : friend Set::Vector operator * (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a, const Set::Matrix3& b);
77 : friend bool operator == (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a, const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& b);
78 : friend Matrix4<AMREX_SPACEDIM, Sym::Diagonal> operator + (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a, const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& b);
79 : friend Matrix4<AMREX_SPACEDIM, Sym::Diagonal> operator - (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a, const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& b);
80 : friend Matrix4<AMREX_SPACEDIM, Sym::Diagonal> operator * (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a, const Set::Scalar& b);
81 : friend Matrix4<AMREX_SPACEDIM, Sym::Diagonal> operator * (const Set::Scalar& b, const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a);
82 : friend Matrix4<AMREX_SPACEDIM, Sym::Diagonal> operator / (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a, const Set::Scalar& b);
83 :
84 :
85 : //AMREX_GPU_HOST_DEVICE void operator = (Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a) {lambda = a.lambda; mu = a.mu;}
86 27593 : AMREX_GPU_HOST_DEVICE void operator += (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a) {
87 82779 : for (int i = 0 ; i < AMREX_SPACEDIM; i++)
88 165558 : for (int j = 0 ; j < AMREX_SPACEDIM; j++)
89 110372 : A[i][j] += a.A[i][j];
90 27593 : }
91 : AMREX_GPU_HOST_DEVICE void operator -= (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a) {
92 : for (int i = 0 ; i < AMREX_SPACEDIM; i++)
93 : for (int j = 0 ; j < AMREX_SPACEDIM; j++)
94 : A[i][j] -= a.A[i][j];
95 : }
96 : AMREX_GPU_HOST_DEVICE void operator *= (const Set::Scalar& alpha) {
97 : for (int i = 0 ; i < AMREX_SPACEDIM; i++)
98 : for (int j = 0 ; j < AMREX_SPACEDIM; j++)
99 : A[i][j] *= alpha;
100 : }
101 : AMREX_GPU_HOST_DEVICE void operator /= (const Set::Scalar& alpha) {
102 : for (int i = 0 ; i < AMREX_SPACEDIM; i++)
103 : for (int j = 0 ; j < AMREX_SPACEDIM; j++)
104 : A[i][j] /= alpha;
105 : }
106 :
107 : Set::Scalar Norm()
108 : {
109 : return getA().lpNorm<2>();
110 : }
111 : bool contains_nan() const
112 : {
113 : if (std::isnan(getA().lpNorm<2>())) return true;
114 : return false;
115 : }
116 :
117 : };
118 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
119 : Set::Matrix operator * (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a, const Set::Matrix& b)
120 : {
121 7682148 : Set::Matrix ret;
122 :
123 23048254 : for (int i = 0; i < AMREX_SPACEDIM; i++)
124 : {
125 46103748 : for (int j = 0; j < AMREX_SPACEDIM; j++)
126 : {
127 30737642 : ret(i, j) = a.A[i][j] * b(i, j);
128 : }
129 : }
130 :
131 7682148 : return ret;
132 : }
133 :
134 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
135 : Set::Vector operator * (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a, const Set::Matrix3& b)
136 : {
137 2519210 : Set::Vector ret = Set::Vector::Zero();
138 :
139 7557630 : for (int i = 0; i < AMREX_SPACEDIM; i++)
140 15115260 : for (int j = 0; j < AMREX_SPACEDIM; j++)
141 20153680 : ret(i) += a.A[i][j] * b(i, j, j);
142 :
143 2519210 : return ret;
144 : }
145 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
146 : bool operator == (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a,
147 : const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& b)
148 : {
149 10 : return a.A == b.A;
150 : }
151 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
152 : Matrix4<AMREX_SPACEDIM, Sym::Diagonal> operator + (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a,
153 : const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& b)
154 : {
155 23424 : Matrix4<AMREX_SPACEDIM, Sym::Diagonal> ret = a;
156 70272 : for (int i = 0; i < AMREX_SPACEDIM; i++)
157 140544 : for (int j = 0; j < AMREX_SPACEDIM; j++)
158 93696 : ret.A[i][j] += b.A[i][j];
159 23424 : return ret;
160 : }
161 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
162 : Matrix4<AMREX_SPACEDIM, Sym::Diagonal> operator - (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a,
163 : const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& b)
164 : {
165 4936104 : Matrix4<AMREX_SPACEDIM, Sym::Diagonal> ret = a;
166 14808312 : for (int i = 0; i < AMREX_SPACEDIM; i++)
167 29616624 : for (int j = 0; j < AMREX_SPACEDIM; j++)
168 19744416 : ret.A[i][j] -= b.A[i][j];
169 4936104 : return ret;
170 : }
171 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
172 : Matrix4<AMREX_SPACEDIM, Sym::Diagonal> operator * (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a,
173 : const Set::Scalar& b)
174 : {
175 4980030 : Matrix4<AMREX_SPACEDIM, Sym::Diagonal> ret = a;
176 14940092 : for (int i = 0; i < AMREX_SPACEDIM; i++)
177 29880192 : for (int j = 0; j < AMREX_SPACEDIM; j++)
178 19920130 : ret.A[i][j] *= b;
179 4980030 : return ret;
180 : }
181 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
182 : Matrix4<AMREX_SPACEDIM, Sym::Diagonal> operator * ( const Set::Scalar& b,
183 : const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a)
184 : {
185 9 : return a*b;
186 : }
187 :
188 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
189 : Matrix4<AMREX_SPACEDIM, Sym::Diagonal> operator / (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a,
190 : const Set::Scalar& b)
191 : {
192 4938792 : Matrix4<AMREX_SPACEDIM, Sym::Diagonal> ret = a;
193 14816376 : for (int i = 0; i < AMREX_SPACEDIM; i++)
194 29632752 : for (int j = 0; j < AMREX_SPACEDIM; j++)
195 19755168 : ret.A[i][j] /= b;
196 4938792 : return ret;
197 : }
198 :
199 :
200 : }
201 : #endif
|
