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 : Set::Matrix getA() const
27 : {
28 : Set::Matrix ret;
29 : for (int i = 0 ; i < AMREX_SPACEDIM; i++)
30 : for (int j = 0 ; j < AMREX_SPACEDIM; j++)
31 : ret(i,j)=A[i][j];
32 : return ret;
33 : }
34 :
35 68811 : 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 : void Print(std::ostream& os)
56 : {
57 : os << "A = \n" << getA();
58 : }
59 8499 : static Matrix4<AMREX_SPACEDIM, Sym::Diagonal> Zero()
60 : {
61 8499 : Matrix4<AMREX_SPACEDIM, Sym::Diagonal> zero;
62 25508 : for (int i = 0 ; i < AMREX_SPACEDIM; i++)
63 51060 : for (int j = 0 ; j < AMREX_SPACEDIM; j++)
64 34051 : zero.A[i][j] = 0.0;
65 8499 : return zero;
66 : }
67 4 : static Matrix4<AMREX_SPACEDIM, Sym::Diagonal> Ones()
68 : {
69 4 : Matrix4<AMREX_SPACEDIM, Sym::Diagonal> ones;
70 12 : for (int i = 0 ; i < AMREX_SPACEDIM; i++)
71 24 : for (int j = 0 ; j < AMREX_SPACEDIM; j++)
72 16 : ones.A[i][j] = 1.0;
73 4 : 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 8463 : AMREX_GPU_HOST_DEVICE void operator += (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a) {
87 25389 : for (int i = 0 ; i < AMREX_SPACEDIM; i++)
88 50778 : for (int j = 0 ; j < AMREX_SPACEDIM; j++)
89 33852 : A[i][j] += a.A[i][j];
90 8463 : }
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 2443539 : Set::Matrix ret;
122 :
123 7332427 : for (int i = 0; i < AMREX_SPACEDIM; i++)
124 : {
125 14672094 : for (int j = 0; j < AMREX_SPACEDIM; j++)
126 : {
127 9783206 : ret(i, j) = a.A[i][j] * b(i, j);
128 : }
129 : }
130 :
131 2443539 : 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 787937 : Set::Vector ret = Set::Vector::Zero();
138 :
139 2363811 : for (int i = 0; i < AMREX_SPACEDIM; i++)
140 4727622 : for (int j = 0; j < AMREX_SPACEDIM; j++)
141 6303496 : ret(i) += a.A[i][j] * b(i, j, j);
142 :
143 787937 : 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 9920 : Matrix4<AMREX_SPACEDIM, Sym::Diagonal> ret = a;
156 29760 : for (int i = 0; i < AMREX_SPACEDIM; i++)
157 59520 : for (int j = 0; j < AMREX_SPACEDIM; j++)
158 39680 : ret.A[i][j] += b.A[i][j];
159 9920 : 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 1541254 : Matrix4<AMREX_SPACEDIM, Sym::Diagonal> ret = a;
166 4623762 : for (int i = 0; i < AMREX_SPACEDIM; i++)
167 9247524 : for (int j = 0; j < AMREX_SPACEDIM; j++)
168 6165016 : ret.A[i][j] -= b.A[i][j];
169 1541254 : 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 1556653 : Matrix4<AMREX_SPACEDIM, Sym::Diagonal> ret = a;
176 4669961 : for (int i = 0; i < AMREX_SPACEDIM; i++)
177 9339930 : for (int j = 0; j < AMREX_SPACEDIM; j++)
178 6226622 : ret.A[i][j] *= b;
179 1556653 : 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 4 : 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 1542390 : Matrix4<AMREX_SPACEDIM, Sym::Diagonal> ret = a;
193 4627170 : for (int i = 0; i < AMREX_SPACEDIM; i++)
194 9254340 : for (int j = 0; j < AMREX_SPACEDIM; j++)
195 6169560 : ret.A[i][j] /= b;
196 1542390 : return ret;
197 : }
198 :
199 :
200 : }
201 : #endif
|
