MFEM  v4.5.2
Finite element discretization library
lor_util.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #ifndef MFEM_LOR_UTIL
13 #define MFEM_LOR_UTIL
14 
15 #include "../../config/config.hpp"
16 #include "../../general/backends.hpp"
17 #include "../../general/globals.hpp"
18 #include "../../linalg/dtensor.hpp"
19 
20 namespace mfem
21 {
22 
23 template <int ORDER>
24 MFEM_HOST_DEVICE inline void LORVertexCoordinates2D(
25  const double *X, int iel_ho, int kx, int ky, double vx[4], double vy[4])
26 {
27  const int dim = 2;
28  const int nd1d = ORDER + 1;
29  const int nvert_per_el = nd1d*nd1d;
30 
31  const int v0 = kx + nd1d*ky;
32  const int v1 = kx + 1 + nd1d*ky;
33  const int v2 = kx + 1 + nd1d*(ky + 1);
34  const int v3 = kx + nd1d*(ky + 1);
35 
36  const int e0 = dim*(v0 + nvert_per_el*iel_ho);
37  const int e1 = dim*(v1 + nvert_per_el*iel_ho);
38  const int e2 = dim*(v2 + nvert_per_el*iel_ho);
39  const int e3 = dim*(v3 + nvert_per_el*iel_ho);
40 
41  // Vertex coordinates
42  vx[0] = X[e0 + 0];
43  vy[0] = X[e0 + 1];
44 
45  vx[1] = X[e1 + 0];
46  vy[1] = X[e1 + 1];
47 
48  vx[2] = X[e2 + 0];
49  vy[2] = X[e2 + 1];
50 
51  vx[3] = X[e3 + 0];
52  vy[3] = X[e3 + 1];
53 }
54 
55 template <int ORDER>
56 MFEM_HOST_DEVICE inline void LORVertexCoordinates3D(
57  const double *X, int iel_ho, int kx, int ky, int kz,
58  double vx[8], double vy[8], double vz[8])
59 {
60  const int dim = 3;
61  const int nd1d = ORDER + 1;
62  const int nvert_per_el = nd1d*nd1d*nd1d;
63 
64  const int v0 = kx + nd1d*(ky + nd1d*kz);
65  const int v1 = kx + 1 + nd1d*(ky + nd1d*kz);
66  const int v2 = kx + 1 + nd1d*(ky + 1 + nd1d*kz);
67  const int v3 = kx + nd1d*(ky + 1 + nd1d*kz);
68  const int v4 = kx + nd1d*(ky + nd1d*(kz + 1));
69  const int v5 = kx + 1 + nd1d*(ky + nd1d*(kz + 1));
70  const int v6 = kx + 1 + nd1d*(ky + 1 + nd1d*(kz + 1));
71  const int v7 = kx + nd1d*(ky + 1 + nd1d*(kz + 1));
72 
73  const int e0 = dim*(v0 + nvert_per_el*iel_ho);
74  const int e1 = dim*(v1 + nvert_per_el*iel_ho);
75  const int e2 = dim*(v2 + nvert_per_el*iel_ho);
76  const int e3 = dim*(v3 + nvert_per_el*iel_ho);
77  const int e4 = dim*(v4 + nvert_per_el*iel_ho);
78  const int e5 = dim*(v5 + nvert_per_el*iel_ho);
79  const int e6 = dim*(v6 + nvert_per_el*iel_ho);
80  const int e7 = dim*(v7 + nvert_per_el*iel_ho);
81 
82  vx[0] = X[e0 + 0];
83  vy[0] = X[e0 + 1];
84  vz[0] = X[e0 + 2];
85 
86  vx[1] = X[e1 + 0];
87  vy[1] = X[e1 + 1];
88  vz[1] = X[e1 + 2];
89 
90  vx[2] = X[e2 + 0];
91  vy[2] = X[e2 + 1];
92  vz[2] = X[e2 + 2];
93 
94  vx[3] = X[e3 + 0];
95  vy[3] = X[e3 + 1];
96  vz[3] = X[e3 + 2];
97 
98  vx[4] = X[e4 + 0];
99  vy[4] = X[e4 + 1];
100  vz[4] = X[e4 + 2];
101 
102  vx[5] = X[e5 + 0];
103  vy[5] = X[e5 + 1];
104  vz[5] = X[e5 + 2];
105 
106  vx[6] = X[e6 + 0];
107  vy[6] = X[e6 + 1];
108  vz[6] = X[e6 + 2];
109 
110  vx[7] = X[e7 + 0];
111  vy[7] = X[e7 + 1];
112  vz[7] = X[e7 + 2];
113 }
114 
115 MFEM_HOST_DEVICE inline void Jacobian2D(
116  const double x, const double y, const double vx[4], const double vy[4],
117  DeviceMatrix &J)
118 {
119  J(0,0) = -(1-y)*vx[0] + (1-y)*vx[1] + y*vx[2] - y*vx[3];
120  J(0,1) = -(1-x)*vx[0] - x*vx[1] + x*vx[2] + (1-x)*vx[3];
121 
122  J(1,0) = -(1-y)*vy[0] + (1-y)*vy[1] + y*vy[2] - y*vy[3];
123  J(1,1) = -(1-x)*vy[0] - x*vy[1] + x*vy[2] + (1-x)*vy[3];
124 }
125 
126 MFEM_HOST_DEVICE inline void Jacobian3D(
127  const double x, const double y, const double z,
128  const double vx[8], const double vy[8], const double vz[8],
129  DeviceMatrix &J)
130 {
131  // c: (1-x)(1-y)(1-z)v0[c] + x (1-y)(1-z)v1[c] + x y (1-z)v2[c] + (1-x) y (1-z)v3[c]
132  // + (1-x)(1-y) z v4[c] + x (1-y) z v5[c] + x y z v6[c] + (1-x) y z v7[c]
133  J(0,0) = -(1-y)*(1-z)*vx[0]
134  + (1-y)*(1-z)*vx[1] + y*(1-z)*vx[2] - y*(1-z)*vx[3]
135  - (1-y)*z*vx[4] + (1-y)*z*vx[5] + y*z*vx[6] - y*z*vx[7];
136 
137  J(0,1) = -(1-x)*(1-z)*vx[0]
138  - x*(1-z)*vx[1] + x*(1-z)*vx[2] + (1-x)*(1-z)*vx[3]
139  - (1-x)*z*vx[4] - x*z*vx[5] + x*z*vx[6] + (1-x)*z*vx[7];
140 
141  J(0,2) = -(1-x)*(1-y)*vx[0] - x*(1-y)*vx[1]
142  - x*y*vx[2] - (1-x)*y*vx[3] + (1-x)*(1-y)*vx[4]
143  + x*(1-y)*vx[5] + x*y*vx[6] + (1-x)*y*vx[7];
144 
145  J(1,0) = -(1-y)*(1-z)*vy[0] + (1-y)*(1-z)*vy[1]
146  + y*(1-z)*vy[2] - y*(1-z)*vy[3] - (1-y)*z*vy[4]
147  + (1-y)*z*vy[5] + y*z*vy[6] - y*z*vy[7];
148 
149  J(1,1) = -(1-x)*(1-z)*vy[0] - x*(1-z)*vy[1]
150  + x*(1-z)*vy[2] + (1-x)*(1-z)*vy[3]- (1-x)*z*vy[4] -
151  x*z*vy[5] + x*z*vy[6] + (1-x)*z*vy[7];
152 
153  J(1,2) = -(1-x)*(1-y)*vy[0] - x*(1-y)*vy[1]
154  - x*y*vy[2] - (1-x)*y*vy[3] + (1-x)*(1-y)*vy[4]
155  + x*(1-y)*vy[5] + x*y*vy[6] + (1-x)*y*vy[7];
156 
157  J(2,0) = -(1-y)*(1-z)*vz[0] + (1-y)*(1-z)*vz[1]
158  + y*(1-z)*vz[2] - y*(1-z)*vz[3]- (1-y)*z*vz[4] +
159  (1-y)*z*vz[5] + y*z*vz[6] - y*z*vz[7];
160 
161  J(2,1) = -(1-x)*(1-z)*vz[0] - x*(1-z)*vz[1]
162  + x*(1-z)*vz[2] + (1-x)*(1-z)*vz[3] - (1-x)*z*vz[4]
163  - x*z*vz[5] + x*z*vz[6] + (1-x)*z*vz[7];
164 
165  J(2,2) = -(1-x)*(1-y)*vz[0] - x*(1-y)*vz[1]
166  - x*y*vz[2] - (1-x)*y*vz[3] + (1-x)*(1-y)*vz[4]
167  + x*(1-y)*vz[5] + x*y*vz[6] + (1-x)*y*vz[7];
168 }
169 
170 MFEM_HOST_DEVICE inline void Adjugate3D(const DeviceMatrix &J, DeviceMatrix &A)
171 {
172  A(0,0) = (J(1,1) * J(2,2)) - (J(1,2) * J(2,1));
173  A(0,1) = (J(2,1) * J(0,2)) - (J(0,1) * J(2,2));
174  A(0,2) = (J(0,1) * J(1,2)) - (J(1,1) * J(0,2));
175  A(1,0) = (J(2,0) * J(1,2)) - (J(1,0) * J(2,2));
176  A(1,1) = (J(0,0) * J(2,2)) - (J(0,2) * J(2,0));
177  A(1,2) = (J(1,0) * J(0,2)) - (J(0,0) * J(1,2));
178  A(2,0) = (J(1,0) * J(2,1)) - (J(2,0) * J(1,1));
179  A(2,1) = (J(2,0) * J(0,1)) - (J(0,0) * J(2,1));
180  A(2,2) = (J(0,0) * J(1,1)) - (J(0,1) * J(1,0));
181 }
182 
183 MFEM_HOST_DEVICE inline double Det2D(DeviceMatrix &J)
184 {
185  return J(0,0)*J(1,1) - J(1,0)*J(0,1);
186 }
187 
188 MFEM_HOST_DEVICE inline double Det3D(DeviceMatrix &J)
189 {
190  return J(0,0) * (J(1,1) * J(2,2) - J(2,1) * J(1,2)) -
191  J(1,0) * (J(0,1) * J(2,2) - J(2,1) * J(0,2)) +
192  J(2,0) * (J(0,1) * J(1,2) - J(1,1) * J(0,2));
193 }
194 
195 }
196 
197 #endif
MFEM_HOST_DEVICE void LORVertexCoordinates3D(const double *X, int iel_ho, int kx, int ky, int kz, double vx[8], double vy[8], double vz[8])
Definition: lor_util.hpp:56
MFEM_HOST_DEVICE double Det3D(DeviceMatrix &J)
Definition: lor_util.hpp:188
MFEM_HOST_DEVICE double Det2D(DeviceMatrix &J)
Definition: lor_util.hpp:183
MFEM_HOST_DEVICE void Jacobian3D(const double x, const double y, const double z, const double vx[8], const double vy[8], const double vz[8], DeviceMatrix &J)
Definition: lor_util.hpp:126
A basic generic Tensor class, appropriate for use on the GPU.
Definition: dtensor.hpp:81
MFEM_HOST_DEVICE void Jacobian2D(const double x, const double y, const double vx[4], const double vy[4], DeviceMatrix &J)
Definition: lor_util.hpp:115
MFEM_HOST_DEVICE void LORVertexCoordinates2D(const double *X, int iel_ho, int kx, int ky, double vx[4], double vy[4])
Definition: lor_util.hpp:24
int dim
Definition: ex24.cpp:53
MFEM_HOST_DEVICE void Adjugate3D(const DeviceMatrix &J, DeviceMatrix &A)
Definition: lor_util.hpp:170