12 #include "../general/forall.hpp" 19 template<
int T_D1D = 0,
int T_Q1D = 0>
20 static void EAMassAssemble1D(
const int NE,
21 const Array<double> &basis,
28 const int D1D = T_D1D ? T_D1D : d1d;
29 const int Q1D = T_Q1D ? T_Q1D : q1d;
30 MFEM_VERIFY(D1D <=
MAX_D1D,
"");
31 MFEM_VERIFY(Q1D <=
MAX_Q1D,
"");
32 auto B =
Reshape(basis.Read(), Q1D, D1D);
33 auto D =
Reshape(padata.Read(), Q1D, NE);
34 auto M =
Reshape(eadata.ReadWrite(), D1D, D1D, NE);
35 MFEM_FORALL_3D(e, NE, D1D, D1D, 1,
37 const int D1D = T_D1D ? T_D1D : d1d;
38 const int Q1D = T_Q1D ? T_Q1D : q1d;
39 constexpr
int MQ1 = T_Q1D ? T_Q1D :
MAX_Q1D;
42 for (
int q = 0; q < Q1D; q++)
44 r_Bi[q] = B(q,MFEM_THREAD_ID(x));
45 r_Bj[q] = B(q,MFEM_THREAD_ID(y));
47 MFEM_FOREACH_THREAD(i1,x,D1D)
49 MFEM_FOREACH_THREAD(j1,y,D1D)
52 for (
int k1 = 0; k1 < Q1D; ++k1)
54 val += r_Bi[k1] * r_Bj[k1] * D(k1, e);
69 template<
int T_D1D = 0,
int T_Q1D = 0>
70 static void EAMassAssemble2D(
const int NE,
71 const Array<double> &basis,
78 const int D1D = T_D1D ? T_D1D : d1d;
79 const int Q1D = T_Q1D ? T_Q1D : q1d;
80 MFEM_VERIFY(D1D <=
MAX_D1D,
"");
81 MFEM_VERIFY(Q1D <=
MAX_Q1D,
"");
82 auto B =
Reshape(basis.Read(), Q1D, D1D);
83 auto D =
Reshape(padata.Read(), Q1D, Q1D, NE);
84 auto M =
Reshape(eadata.ReadWrite(), D1D, D1D, D1D, D1D, NE);
85 MFEM_FORALL_3D(e, NE, D1D, D1D, 1,
87 const int D1D = T_D1D ? T_D1D : d1d;
88 const int Q1D = T_Q1D ? T_Q1D : q1d;
89 constexpr
int MD1 = T_D1D ? T_D1D :
MAX_D1D;
90 constexpr
int MQ1 = T_Q1D ? T_Q1D :
MAX_Q1D;
92 for (
int d = 0; d < D1D; d++)
94 for (
int q = 0; q < Q1D; q++)
99 MFEM_SHARED
double s_D[MQ1][MQ1];
100 MFEM_FOREACH_THREAD(k1,x,Q1D)
102 MFEM_FOREACH_THREAD(k2,y,Q1D)
104 s_D[k1][k2] = D(k1,k2,e);
108 MFEM_FOREACH_THREAD(i1,x,D1D)
110 MFEM_FOREACH_THREAD(i2,y,D1D)
112 for (
int j1 = 0; j1 < D1D; ++j1)
114 for (
int j2 = 0; j2 < D1D; ++j2)
117 for (
int k1 = 0; k1 < Q1D; ++k1)
119 for (
int k2 = 0; k2 < Q1D; ++k2)
121 val += r_B[k1][i1] * r_B[k1][j1]
122 * r_B[k2][i2] * r_B[k2][j2]
128 M(i1, i2, j1, j2, e) += val;
132 M(i1, i2, j1, j2, e) = val;
141 template<
int T_D1D = 0,
int T_Q1D = 0>
142 static void EAMassAssemble3D(
const int NE,
143 const Array<double> &basis,
144 const Vector &padata,
150 const int D1D = T_D1D ? T_D1D : d1d;
151 const int Q1D = T_Q1D ? T_Q1D : q1d;
152 MFEM_VERIFY(D1D <=
MAX_D1D,
"");
153 MFEM_VERIFY(Q1D <=
MAX_Q1D,
"");
154 auto B =
Reshape(basis.Read(), Q1D, D1D);
155 auto D =
Reshape(padata.Read(), Q1D, Q1D, Q1D, NE);
156 auto M =
Reshape(eadata.ReadWrite(), D1D, D1D, D1D, D1D, D1D, D1D, NE);
157 MFEM_FORALL_3D(e, NE, D1D, D1D, D1D,
159 const int D1D = T_D1D ? T_D1D : d1d;
160 const int Q1D = T_Q1D ? T_Q1D : q1d;
161 constexpr
int MD1 = T_D1D ? T_D1D :
MAX_D1D;
162 constexpr
int MQ1 = T_Q1D ? T_Q1D :
MAX_Q1D;
163 double r_B[MQ1][MD1];
164 for (
int d = 0; d < D1D; d++)
166 for (
int q = 0; q < Q1D; q++)
171 MFEM_SHARED
double s_D[MQ1][MQ1][MQ1];
172 MFEM_FOREACH_THREAD(k1,x,Q1D)
174 MFEM_FOREACH_THREAD(k2,y,Q1D)
176 MFEM_FOREACH_THREAD(k3,z,Q1D)
178 s_D[k1][k2][k3] = D(k1,k2,k3,e);
183 MFEM_FOREACH_THREAD(i1,x,D1D)
185 MFEM_FOREACH_THREAD(i2,y,D1D)
187 MFEM_FOREACH_THREAD(i3,z,D1D)
189 for (
int j1 = 0; j1 < D1D; ++j1)
191 for (
int j2 = 0; j2 < D1D; ++j2)
193 for (
int j3 = 0; j3 < D1D; ++j3)
196 for (
int k1 = 0; k1 < Q1D; ++k1)
198 for (
int k2 = 0; k2 < Q1D; ++k2)
200 for (
int k3 = 0; k3 < Q1D; ++k3)
202 val += r_B[k1][i1] * r_B[k1][j1]
203 * r_B[k2][i2] * r_B[k2][j2]
204 * r_B[k3][i3] * r_B[k3][j3]
211 M(i1, i2, i3, j1, j2, j3, e) += val;
215 M(i1, i2, i3, j1, j2, j3, e) = val;
237 case 0x22:
return EAMassAssemble1D<2,2>(
ne,B,
pa_data,ea_data,
add);
238 case 0x33:
return EAMassAssemble1D<3,3>(
ne,B,
pa_data,ea_data,
add);
239 case 0x44:
return EAMassAssemble1D<4,4>(
ne,B,
pa_data,ea_data,
add);
240 case 0x55:
return EAMassAssemble1D<5,5>(
ne,B,
pa_data,ea_data,
add);
241 case 0x66:
return EAMassAssemble1D<6,6>(
ne,B,
pa_data,ea_data,
add);
242 case 0x77:
return EAMassAssemble1D<7,7>(
ne,B,
pa_data,ea_data,
add);
243 case 0x88:
return EAMassAssemble1D<8,8>(
ne,B,
pa_data,ea_data,
add);
244 case 0x99:
return EAMassAssemble1D<9,9>(
ne,B,
pa_data,ea_data,
add);
245 default:
return EAMassAssemble1D(
ne,B,
pa_data,ea_data,
add,
253 case 0x22:
return EAMassAssemble2D<2,2>(
ne,B,
pa_data,ea_data,
add);
254 case 0x33:
return EAMassAssemble2D<3,3>(
ne,B,
pa_data,ea_data,
add);
255 case 0x44:
return EAMassAssemble2D<4,4>(
ne,B,
pa_data,ea_data,
add);
256 case 0x55:
return EAMassAssemble2D<5,5>(
ne,B,
pa_data,ea_data,
add);
257 case 0x66:
return EAMassAssemble2D<6,6>(
ne,B,
pa_data,ea_data,
add);
258 case 0x77:
return EAMassAssemble2D<7,7>(
ne,B,
pa_data,ea_data,
add);
259 case 0x88:
return EAMassAssemble2D<8,8>(
ne,B,
pa_data,ea_data,
add);
260 case 0x99:
return EAMassAssemble2D<9,9>(
ne,B,
pa_data,ea_data,
add);
261 default:
return EAMassAssemble2D(
ne,B,
pa_data,ea_data,
add,
269 case 0x23:
return EAMassAssemble3D<2,3>(
ne,B,
pa_data,ea_data,
add);
270 case 0x34:
return EAMassAssemble3D<3,4>(
ne,B,
pa_data,ea_data,
add);
271 case 0x45:
return EAMassAssemble3D<4,5>(
ne,B,
pa_data,ea_data,
add);
272 case 0x56:
return EAMassAssemble3D<5,6>(
ne,B,
pa_data,ea_data,
add);
273 case 0x67:
return EAMassAssemble3D<6,7>(
ne,B,
pa_data,ea_data,
add);
274 case 0x78:
return EAMassAssemble3D<7,8>(
ne,B,
pa_data,ea_data,
add);
275 case 0x89:
return EAMassAssemble3D<8,9>(
ne,B,
pa_data,ea_data,
add);
276 default:
return EAMassAssemble3D(
ne,B,
pa_data,ea_data,
add,
280 MFEM_ABORT(
"Unknown kernel.");
virtual void AssemblePA(const FiniteElementSpace &fes)
Method defining partial assembly.
void add(const Vector &v1, const Vector &v2, Vector &v)
const DofToQuad * maps
Not owned.
virtual void AssembleEA(const FiniteElementSpace &fes, Vector &emat, const bool add)
Method defining element assembly.
Mesh * GetMesh() const
Returns the mesh.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Array< double > B
Basis functions evaluated at quadrature points.
int GetNE() const
Returns number of elements.
MFEM_HOST_DEVICE DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims... dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.