omega_h
Reliable mesh adaptation
Omega_h_matrix.hpp
1 #ifndef OMEGA_H_MATRIX_HPP
2 #define OMEGA_H_MATRIX_HPP
3 
4 #include <Omega_h_array.hpp>
5 #include <Omega_h_vector.hpp>
6 
7 namespace Omega_h {
8 
9 /* column-first storage and indexing !!! */
10 
11 #ifdef OMEGA_H_USE_KOKKOS
12 
13 template <Int m, Int n>
14 class Matrix : public Few<Vector<m>, n> {
15  public:
16  OMEGA_H_INLINE Matrix() {}
17  /* these constructors accept the matrix in
18  * row-first order for convenience.
19  */
20  inline Matrix(std::initializer_list<Vector<m>> l) : Few<Vector<m>, n>(l) {}
21  inline Matrix(std::initializer_list<Real> l);
22  OMEGA_H_INLINE Matrix(Matrix<m, n> const& rhs) : Few<Vector<m>, n>(rhs) {}
23  OMEGA_H_INLINE Real& operator()(Int i, Int j) {
24  return Few<Vector<m>, n>::operator[](j)[i];
25  }
26  OMEGA_H_INLINE Real const& operator()(Int i, Int j) const {
27  return Few<Vector<m>, n>::operator[](j)[i];
28  }
29 };
30 
31 #else
32 
33 template <Int m, Int n>
34 class Matrix : public Few<Vector<m>, n> {
35  public:
36  inline Matrix() = default;
37  /* these constructors accept the matrix in
38  * row-first order for convenience.
39  */
40  inline Matrix(std::initializer_list<Vector<m>> l) : Few<Vector<m>, n>(l) {}
41  inline Matrix(std::initializer_list<Real> l);
42  inline Matrix(Matrix const&) = default;
43  inline Matrix(Matrix&&) = default;
44  inline Matrix& operator=(Matrix const&) = default;
45  inline Matrix& operator=(Matrix&&) = default;
46  OMEGA_H_INLINE Real& operator()(Int i, Int j) OMEGA_H_NOEXCEPT {
47  return Few<Vector<m>, n>::operator[](j)[i];
48  }
49  OMEGA_H_INLINE Real const& operator()(Int i, Int j) const OMEGA_H_NOEXCEPT {
50  return Few<Vector<m>, n>::operator[](j)[i];
51  }
52 };
53 
54 #endif
55 
56 template <Int dim>
57 using Tensor = Matrix<dim, dim>;
58 
59 template <Int m, Int n>
60 OMEGA_H_INLINE Real* scalar_ptr(Matrix<m, n>& a) {
61  return &a[0][0];
62 }
63 template <Int m, Int n>
64 OMEGA_H_INLINE Real const* scalar_ptr(Matrix<m, n> const& a) {
65  return &a[0][0];
66 }
67 
68 template <Int m, Int n>
69 inline Matrix<m, n>::Matrix(std::initializer_list<Real> l) {
70  Int k = 0;
71  for (Real v : l) {
72  (*this)[k % n][k / n] = v;
73  ++k;
74  }
75 }
76 
77 template <Int m, Int n>
78 OMEGA_H_INLINE Vector<m> operator*(
79  Matrix<m, n> a, Vector<n> b)OMEGA_H_NOEXCEPT {
80  Vector<m> c = a[0] * b[0];
81  for (Int j = 1; j < n; ++j) c = c + a[j] * b[j];
82  return c;
83 }
84 
85 template <Int m, Int n, Int p>
86 OMEGA_H_INLINE Matrix<m, n> operator*(Matrix<m, p> a, Matrix<p, n> b) {
87  Matrix<m, n> c;
88  for (Int j = 0; j < n; ++j) c[j] = a * b[j];
89  return c;
90 }
91 
92 template <Int m, Int n>
93 OMEGA_H_INLINE Matrix<n, m> transpose(Matrix<m, n> a) {
94  Matrix<n, m> b;
95  for (Int i = 0; i < m; ++i)
96  for (Int j = 0; j < n; ++j) b[i][j] = a[j][i];
97  return b;
98 }
99 
100 template <Int max_m, Int max_n>
101 OMEGA_H_INLINE Matrix<max_m, max_n> identity_matrix(Int m, Int n) {
102  Matrix<max_m, max_n> a;
103  for (Int j = 0; j < n; ++j)
104  for (Int i = 0; i < m; ++i) a[j][i] = (i == j);
105  return a;
106 }
107 
108 template <Int max_m, Int max_n>
109 OMEGA_H_INLINE Matrix<max_m, max_n> identity_matrix() {
110  return identity_matrix<max_m, max_n>(max_m, max_n);
111 }
112 
113 template <Int dim>
114 OMEGA_H_INLINE Tensor<dim> identity_tensor() {
115  return identity_matrix<dim, dim>();
116 }
117 
118 OMEGA_H_INLINE Matrix<1, 1> matrix_1x1(Real a) {
119  Matrix<1, 1> o;
120  o[0][0] = a;
121  return o;
122 }
123 
124 OMEGA_H_INLINE Tensor<1> tensor_1(Real a) { return matrix_1x1(a); }
125 
126 OMEGA_H_INLINE Matrix<2, 2> matrix_2x2(Real a, Real b, Real c, Real d) {
127  Matrix<2, 2> o;
128  o[0] = vector_2(a, c);
129  o[1] = vector_2(b, d);
130  return o;
131 }
132 
133 OMEGA_H_INLINE Tensor<2> tensor_2(Real a, Real b, Real c, Real d) {
134  return matrix_2x2(a, b, c, d);
135 }
136 
137 OMEGA_H_INLINE Matrix<3, 3> matrix_3x3(
138  Real a, Real b, Real c, Real d, Real e, Real f, Real g, Real h, Real i) {
139  Matrix<3, 3> o;
140  o[0] = vector_3(a, d, g);
141  o[1] = vector_3(b, e, h);
142  o[2] = vector_3(c, f, i);
143  return o;
144 }
145 
146 OMEGA_H_INLINE Tensor<3> tensor_3(
147  Real a, Real b, Real c, Real d, Real e, Real f, Real g, Real h, Real i) {
148  return matrix_3x3(a, b, c, d, e, f, g, h, i);
149 }
150 
151 template <Int m, Int n>
152 OMEGA_H_INLINE Matrix<m, n> operator*(Matrix<m, n> a, Real b) {
153  Matrix<m, n> c;
154  for (Int j = 0; j < n; ++j) c[j] = a[j] * b;
155  return c;
156 }
157 
158 template <Int m, Int n>
159 OMEGA_H_INLINE Matrix<m, n>& operator*=(Matrix<m, n>& a, Real b) {
160  a = a * b;
161  return a;
162 }
163 
164 template <Int m, Int n>
165 OMEGA_H_INLINE Matrix<m, n> operator*(Real a, Matrix<m, n> b) {
166  return b * a;
167 }
168 
169 template <Int m, Int n>
170 OMEGA_H_INLINE Matrix<m, n> operator/(Matrix<m, n> a, Real b) {
171  Matrix<m, n> c;
172  for (Int j = 0; j < n; ++j) c[j] = a[j] / b;
173  return c;
174 }
175 
176 template <Int m, Int n>
177 OMEGA_H_INLINE Matrix<m, n>& operator/=(Matrix<m, n>& a, Real b) {
178  a = a / b;
179  return a;
180 }
181 
182 template <Int m, Int n>
183 OMEGA_H_INLINE Matrix<m, n> operator+(Matrix<m, n> a, Matrix<m, n> b) {
184  Matrix<m, n> c;
185  for (Int j = 0; j < n; ++j) c[j] = a[j] + b[j];
186  return c;
187 }
188 
189 template <Int m, Int n>
190 OMEGA_H_INLINE Matrix<m, n>& operator+=(Matrix<m, n>& a, Matrix<m, n> b) {
191  a = a + b;
192  return a;
193 }
194 
195 template <Int m, Int n>
196 OMEGA_H_INLINE Matrix<m, n> operator-(Matrix<m, n> a, Matrix<m, n> b) {
197  Matrix<m, n> c;
198  for (Int j = 0; j < n; ++j) c[j] = a[j] - b[j];
199  return c;
200 }
201 
202 template <Int m, Int n>
203 OMEGA_H_INLINE Matrix<m, n>& operator-=(Matrix<m, n>& a, Matrix<m, n> b) {
204  a = a - b;
205  return a;
206 }
207 
208 template <Int m, Int n>
209 OMEGA_H_INLINE Matrix<m, n> operator-(Matrix<m, n> a) {
210  Matrix<m, n> c;
211  for (Int j = 0; j < n; ++j) c[j] = -a[j];
212  return c;
213 }
214 
215 template <Int m, Int n>
216 OMEGA_H_INLINE Real max_norm(Matrix<m, n> a) {
217  Real x = 0.0;
218  for (Int j = 0; j < n; ++j)
219  for (Int i = 0; i < m; ++i) x = max2(x, std::abs(a[j][i]));
220  return x;
221 }
222 
223 template <Int m, Int n>
224 OMEGA_H_INLINE bool are_close(
225  Matrix<m, n> a, Matrix<m, n> b, Real tol = EPSILON, Real floor = EPSILON) {
226  for (Int j = 0; j < n; ++j)
227  if (!are_close(a[j], b[j], tol, floor)) return false;
228  return true;
229 }
230 
231 template <Int m, Int n>
232 OMEGA_H_INLINE Matrix<m, n> outer_product(Vector<m> a, Vector<n> b) {
233  Matrix<m, n> c;
234  for (Int j = 0; j < n; ++j)
235  for (Int i = 0; i < m; ++i) c[j][i] = a[i] * b[j];
236  return c;
237 }
238 
239 template <Int m>
240 OMEGA_H_INLINE Real trace(Tensor<m> a) {
241  Real t = a[0][0];
242  for (Int i = 1; i < m; ++i) t += a[i][i];
243  return t;
244 }
245 
246 template <Int m>
247 OMEGA_H_INLINE Vector<m> diagonal(Tensor<m> a) {
248  Vector<m> v;
249  for (Int i = 0; i < m; ++i) v[i] = a[i][i];
250  return v;
251 }
252 
253 template <Int m, Int n>
254 OMEGA_H_INLINE Matrix<m, n> zero_matrix() {
255  Matrix<m, n> a;
256  for (Int j = 0; j < n; ++j) a[j] = zero_vector<m>();
257  return a;
258 }
259 
260 OMEGA_H_INLINE Real determinant(Tensor<1> m) { return m[0][0]; }
261 
262 OMEGA_H_INLINE Real determinant(Tensor<2> m) {
263  Real a = m[0][0];
264  Real b = m[1][0];
265  Real c = m[0][1];
266  Real d = m[1][1];
267  return a * d - b * c;
268 }
269 
270 OMEGA_H_INLINE Real determinant(Tensor<3> m) {
271  Real a = m[0][0];
272  Real b = m[1][0];
273  Real c = m[2][0];
274  Real d = m[0][1];
275  Real e = m[1][1];
276  Real f = m[2][1];
277  Real g = m[0][2];
278  Real h = m[1][2];
279  Real i = m[2][2];
280  return (a * e * i) + (b * f * g) + (c * d * h) - (c * e * g) - (b * d * i) -
281  (a * f * h);
282 }
283 
284 template <Int max_m, Int max_n>
285 OMEGA_H_INLINE Real inner_product(
286  Int m, Int n, Matrix<max_m, max_n> a, Matrix<max_m, max_n> b) {
287  Real out = 0.0;
288  for (Int j = 0; j < n; ++j) {
289  for (Int i = 0; i < m; ++i) {
290  out += a[j][i] * b[j][i];
291  }
292  }
293  return out;
294 }
295 
296 template <Int m, Int n>
297 OMEGA_H_INLINE Real inner_product(Matrix<m, n> a, Matrix<m, n> b) {
298  return inner_product(m, n, a, b);
299 }
300 
301 template <Int max_m, Int max_n>
302 OMEGA_H_INLINE Real norm(
303  Int const m, Int const n, Matrix<max_m, max_n> const & a) {
304  return std::sqrt(inner_product(m, n, a, a));
305 }
306 
307 template <Int m, Int n>
308 OMEGA_H_INLINE Real norm(Matrix<m, n> const & a) {
309  return norm(m, n, a);
310 }
311 
312 OMEGA_H_INLINE Tensor<3> cross(Vector<3> const & a) {
313  return matrix_3x3(0, -a[2], a[1], a[2], 0, -a[0], -a[1], a[0], 0);
314 }
315 
316 OMEGA_H_INLINE Vector<3> uncross(Tensor<3> const & c) {
317  return 0.5 *
318  vector_3(c[1][2] - c[2][1], c[2][0] - c[0][2], c[0][1] - c[1][0]);
319 }
320 
321 OMEGA_H_INLINE Tensor<1> invert(Tensor<1> const & m) {
322  return matrix_1x1(1.0 / m[0][0]);
323 }
324 
325 OMEGA_H_INLINE Tensor<2> invert(Tensor<2> const & m) {
326  Real a = m[0][0];
327  Real b = m[1][0];
328  Real c = m[0][1];
329  Real d = m[1][1];
330  return matrix_2x2(d, -b, -c, a) / determinant(m);
331 }
332 
333 OMEGA_H_INLINE Tensor<3> invert(Tensor<3> const & a) {
334  Tensor<3> b;
335  b[0] = cross(a[1], a[2]);
336  b[1] = cross(a[2], a[0]);
337  b[2] = cross(a[0], a[1]);
338  return transpose(b) / determinant(a);
339 }
340 
341 template <Int m, Int n>
342 OMEGA_H_INLINE typename std::enable_if<(n == m), Matrix<n, m>>::type
343 pseudo_invert(Matrix<m, n> const & a) {
344  return invert(a);
345 }
346 
347 template <Int m, Int n>
348 OMEGA_H_INLINE typename std::enable_if<(n < m), Matrix<n, m>>::type
349 pseudo_invert(Matrix<m, n> const &a) {
350  auto at = transpose(a);
351  return invert(at * a) * at;
352 }
353 
354 template <Int m, Int n>
355 OMEGA_H_INLINE typename std::enable_if<(n > m), Matrix<n, m>>::type
356 pseudo_invert(Matrix<m, n> const &a) {
357  auto at = transpose(a);
358  return at * invert(a * at);
359 }
360 
361 template <Int m>
362 OMEGA_H_INLINE Tensor<m> diagonal(Vector<m> v) {
363  Tensor<m> a;
364  for (Int i = 0; i < m; ++i)
365  for (Int j = i + 1; j < m; ++j) a[i][j] = a[j][i] = 0.0;
366  for (Int i = 0; i < m; ++i) a[i][i] = v[i];
367  return a;
368 }
369 
370 OMEGA_H_INLINE constexpr Int symm_ncomps(Int dim) {
371  return ((dim + 1) * dim) / 2;
372 }
373 
374 OMEGA_H_INLINE Vector<1> symm2vector(Tensor<1> symm) {
375  return vector_1(symm[0][0]);
376 }
377 
378 OMEGA_H_INLINE Tensor<1> vector2symm(Vector<1> v) { return matrix_1x1(v[0]); }
379 
380 OMEGA_H_INLINE Vector<3> symm2vector(Tensor<2> symm) {
381  Vector<3> v;
382  v[0] = symm[0][0];
383  v[1] = symm[1][1];
384  v[2] = symm[1][0];
385  return v;
386 }
387 
388 OMEGA_H_INLINE Tensor<2> vector2symm(Vector<3> v) {
389  Tensor<2> symm;
390  symm[0][0] = v[0];
391  symm[1][1] = v[1];
392  symm[1][0] = v[2];
393  symm[0][1] = symm[1][0];
394  return symm;
395 }
396 
397 OMEGA_H_INLINE Vector<6> symm2vector(Tensor<3> symm) {
398  Vector<6> v;
399  v[0] = symm[0][0];
400  v[1] = symm[1][1];
401  v[2] = symm[2][2];
402  v[3] = symm[1][0];
403  v[4] = symm[2][1];
404  v[5] = symm[2][0];
405  return v;
406 }
407 
408 OMEGA_H_INLINE Tensor<3> vector2symm(Vector<6> v) {
409  Tensor<3> symm;
410  symm[0][0] = v[0];
411  symm[1][1] = v[1];
412  symm[2][2] = v[2];
413  symm[1][0] = v[3];
414  symm[2][1] = v[4];
415  symm[2][0] = v[5];
416  symm[0][1] = symm[1][0];
417  symm[1][2] = symm[2][1];
418  symm[0][2] = symm[2][0];
419  return symm;
420 }
421 
422 /* Symmetric metric tensor storage convention used by
423  INRIA: https://hal.inria.fr/inria-00363007/document */
424 OMEGA_H_INLINE Tensor<1> vector2symm_inria(Vector<1> v) {
425  return matrix_1x1(v[0]);
426 }
427 
428 OMEGA_H_INLINE Tensor<2> vector2symm_inria(Vector<3> v) {
429  Tensor<2> symm;
430  symm[0][0] = v[0];
431  symm[0][1] = v[1];
432  symm[1][1] = v[2];
433  symm[1][0] = symm[0][1];
434  return symm;
435 }
436 
437 OMEGA_H_INLINE Tensor<3> vector2symm_inria(Vector<6> v) {
438  Tensor<3> symm;
439  symm[0][0] = v[0];
440  symm[0][1] = v[1];
441  symm[1][1] = v[2];
442  symm[0][2] = v[3];
443  symm[1][2] = v[4];
444  symm[2][2] = v[5];
445  symm[1][0] = symm[0][1];
446  symm[2][0] = symm[0][2];
447  symm[2][1] = symm[1][2];
448  return symm;
449 }
450 
451 OMEGA_H_INLINE Vector<1> symm2vector_inria(Tensor<1> symm) {
452  return vector_1(symm[0][0]);
453 }
454 
455 OMEGA_H_INLINE Vector<3> symm2vector_inria(Tensor<2> symm) {
456  Vector<3> v;
457  v[0] = symm[0][0];
458  v[1] = symm[0][1];
459  v[2] = symm[1][1];
460  return v;
461 }
462 
463 OMEGA_H_INLINE Vector<6> symm2vector_inria(Tensor<3> symm) {
464  Vector<6> v;
465  v[0] = symm[0][0];
466  v[1] = symm[0][1];
467  v[2] = symm[1][1];
468  v[3] = symm[0][2];
469  v[4] = symm[1][2];
470  v[5] = symm[2][2];
471  return v;
472 }
473 
474 OMEGA_H_INLINE constexpr Int matrix_ncomps(Int m, Int n) { return m * n; }
475 
476 template <Int m, Int n>
477 OMEGA_H_INLINE Vector<matrix_ncomps(m, n)> matrix2vector(Matrix<m, n> a) {
478  Vector<matrix_ncomps(m, n)> v;
479  for (Int i = 0; i < m; ++i) {
480  for (Int j = 0; j < n; ++j) {
481  v[i * n + j] = a(i, j);
482  }
483  }
484  return v;
485 }
486 
487 template <Int m, Int n>
488 OMEGA_H_INLINE Matrix<m, n> vector2matrix(Vector<matrix_ncomps(m, n)> v) {
489  Matrix<m, n> a;
490  for (Int i = 0; i < m; ++i) {
491  for (Int j = 0; j < n; ++j) {
492  a(i, j) = v[i * n + j];
493  }
494  }
495  return a;
496 }
497 
498 template <Int n>
499 OMEGA_H_DEVICE void set_symm(Write<Real> const& a, Int i, Tensor<n> symm) {
500  set_vector(a, i, symm2vector(symm));
501 }
502 
503 template <Int n, typename Arr>
504 OMEGA_H_DEVICE Tensor<n> get_symm(Arr const& a, Int i) {
505  return vector2symm(get_vector<symm_ncomps(n)>(a, i));
506 }
507 
508 template <Int m, Int n>
509 OMEGA_H_DEVICE void set_matrix(
510  Write<Real> const& matrices, Int i, Matrix<m, n> matrix) {
511  set_vector(matrices, i, matrix2vector(matrix));
512 }
513 
514 template <Int m, Int n, typename Arr>
515 OMEGA_H_DEVICE Matrix<m, n> get_matrix(Arr const& matrices, Int const i) {
516  return vector2matrix<m, n>(get_vector<matrix_ncomps(m, n)>(matrices, i));
517 }
518 
519 /* Rodrigues' Rotation Formula */
520 OMEGA_H_INLINE Tensor<3> rotate(Real const angle, Vector<3> const axis) {
521  return std::cos(angle) * identity_matrix<3, 3>() +
522  std::sin(angle) * cross(axis) +
523  (1 - std::cos(angle)) * outer_product(axis, axis);
524 }
525 
526 OMEGA_H_INLINE Tensor<2> rotate(Real const angle) {
527  return matrix_2x2(
528  std::cos(angle), -std::sin(angle), std::sin(angle), std::cos(angle));
529 }
530 
531 OMEGA_H_INLINE Real rotation_angle(Tensor<2> const r) {
532  auto const cos_theta = 0.5 * trace(r);
533  auto const sin_theta = 0.5 * (r(1, 0) - r(0, 1));
534  return std::atan2(sin_theta, cos_theta);
535 }
536 
537 OMEGA_H_INLINE Real rotation_angle(Tensor<3> const r) {
538  auto const cos_theta = 0.5 * (trace(r) - 1.0);
539  return std::acos(cos_theta);
540 }
541 
542 OMEGA_H_INLINE Tensor<1> form_ortho_basis(Vector<1> const v) {
543  return tensor_1(v[0]);
544 }
545 
546 OMEGA_H_INLINE Tensor<2> form_ortho_basis(Vector<2> const v) {
547  Tensor<2> A;
548  A[0] = v;
549  A[1] = perp(A[0]);
550  return A;
551 }
552 
553 /* Duff, Tom, et al.
554  "Building an orthonormal basis, revisited."
555  Journal of Computer Graphics Techniques Vol 6.1 (2017). */
556 OMEGA_H_INLINE Tensor<3> form_ortho_basis(Vector<3> const v) {
557  Tensor<3> A;
558  A[0] = v;
559  auto sign = std::copysign(1.0, v(2));
560  const auto a = -1.0 / (sign + v(2));
561  const auto b = v(0) * v(1) * a;
562  A[1] = vector_3(1.0 + sign * v(0) * v(0) * a, sign * b, -sign * v(0));
563  A[2] = vector_3(b, sign + v(1) * v(1) * a, -v(1));
564  return A;
565 }
566 
567 template <Int dim>
568 OMEGA_H_INLINE Tensor<dim> deviator(Tensor<dim> const a) {
569  return (a - ((1.0 / dim) * trace(a) * identity_matrix<dim, dim>()));
570 }
571 
572 template <int new_dim, int old_dim>
573 OMEGA_H_INLINE Tensor<new_dim> resize(Tensor<old_dim> const m) {
574  constexpr int min_dim = Omega_h::min2(new_dim, old_dim);
575  Tensor<new_dim> m2;
576  for (int i = 0; i < min_dim; ++i)
577  for (int j = 0; j < min_dim; ++j) m2(i, j) = m(i, j);
578  for (int i = 0; i < new_dim; ++i)
579  for (int j = min_dim; j < new_dim; ++j) m2(i, j) = m2(j, i) = 0.0;
580  return m2;
581 }
582 
583 template <Int dim>
584 Reals repeat_symm(LO const n, Tensor<dim> const symm);
585 extern template Reals repeat_symm(LO const n, Tensor<3> const symm);
586 extern template Reals repeat_symm(LO const n, Tensor<2> const symm);
587 extern template Reals repeat_symm(LO const n, Tensor<1> const symm);
588 
589 Reals resize_symms(Reals old_symms, Int old_dim, Int new_dim);
590 
591 template <Int dim>
592 Reals repeat_matrix(LO const n, Tensor<dim> m);
593 
594 extern template Reals repeat_matrix(LO const n, Tensor<3> const m);
595 extern template Reals repeat_matrix(LO const n, Tensor<2> const m);
596 extern template Reals repeat_matrix(LO const n, Tensor<1> const m);
597 
598 Reals matrices_times_vectors(Reals ms, Reals vs, Int dim);
599 Reals matrices_times_matrices(Reals ms, Reals vs, Int dim);
600 
601 Reals symms_inria2osh(Int dim, Reals symms);
602 Reals symms_osh2inria(Int dim, Reals symms);
603 
604 Reals matrices_to_symms(Reals const matrices, Int const dim);
605 
606 } // end namespace Omega_h
607 
608 #endif
Omega_h::Few
Definition: Omega_h_few.hpp:61
Omega_h::Matrix
Definition: Omega_h_matrix.hpp:34
Omega_h
Definition: amr_mpi_test.cpp:6
Generated by doxygen 1.9.1