MFEM  v4.3.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
bilininteg_mass_ea.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2021, 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 #include "../general/forall.hpp"
13 #include "bilininteg.hpp"
14 #include "gridfunc.hpp"
15 
16 namespace mfem
17 {
18 
19 template<int T_D1D = 0, int T_Q1D = 0>
20 static void EAMassAssemble1D(const int NE,
21  const Array<double> &basis,
22  const Vector &padata,
23  Vector &eadata,
24  const bool add,
25  const int d1d = 0,
26  const int q1d = 0)
27 {
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,
36  {
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;
40  double r_Bi[MQ1];
41  double r_Bj[MQ1];
42  for (int q = 0; q < Q1D; q++)
43  {
44  r_Bi[q] = B(q,MFEM_THREAD_ID(x));
45  r_Bj[q] = B(q,MFEM_THREAD_ID(y));
46  }
47  MFEM_FOREACH_THREAD(i1,x,D1D)
48  {
49  MFEM_FOREACH_THREAD(j1,y,D1D)
50  {
51  double val = 0.0;
52  for (int k1 = 0; k1 < Q1D; ++k1)
53  {
54  val += r_Bi[k1] * r_Bj[k1] * D(k1, e);
55  }
56  if (add)
57  {
58  M(i1, j1, e) += val;
59  }
60  else
61  {
62  M(i1, j1, e) = val;
63  }
64  }
65  }
66  });
67 }
68 
69 template<int T_D1D = 0, int T_Q1D = 0>
70 static void EAMassAssemble2D(const int NE,
71  const Array<double> &basis,
72  const Vector &padata,
73  Vector &eadata,
74  const bool add,
75  const int d1d = 0,
76  const int q1d = 0)
77 {
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,
86  {
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;
91  double r_B[MQ1][MD1];
92  for (int d = 0; d < D1D; d++)
93  {
94  for (int q = 0; q < Q1D; q++)
95  {
96  r_B[q][d] = B(q,d);
97  }
98  }
99  MFEM_SHARED double s_D[MQ1][MQ1];
100  MFEM_FOREACH_THREAD(k1,x,Q1D)
101  {
102  MFEM_FOREACH_THREAD(k2,y,Q1D)
103  {
104  s_D[k1][k2] = D(k1,k2,e);
105  }
106  }
107  MFEM_SYNC_THREAD;
108  MFEM_FOREACH_THREAD(i1,x,D1D)
109  {
110  MFEM_FOREACH_THREAD(i2,y,D1D)
111  {
112  for (int j1 = 0; j1 < D1D; ++j1)
113  {
114  for (int j2 = 0; j2 < D1D; ++j2)
115  {
116  double val = 0.0;
117  for (int k1 = 0; k1 < Q1D; ++k1)
118  {
119  for (int k2 = 0; k2 < Q1D; ++k2)
120  {
121  val += r_B[k1][i1] * r_B[k1][j1]
122  * r_B[k2][i2] * r_B[k2][j2]
123  * s_D[k1][k2];
124  }
125  }
126  if (add)
127  {
128  M(i1, i2, j1, j2, e) += val;
129  }
130  else
131  {
132  M(i1, i2, j1, j2, e) = val;
133  }
134  }
135  }
136  }
137  }
138  });
139 }
140 
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,
145  Vector &eadata,
146  const bool add,
147  const int d1d = 0,
148  const int q1d = 0)
149 {
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,
158  {
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++)
165  {
166  for (int q = 0; q < Q1D; q++)
167  {
168  r_B[q][d] = B(q,d);
169  }
170  }
171  MFEM_SHARED double s_D[MQ1][MQ1][MQ1];
172  MFEM_FOREACH_THREAD(k1,x,Q1D)
173  {
174  MFEM_FOREACH_THREAD(k2,y,Q1D)
175  {
176  MFEM_FOREACH_THREAD(k3,z,Q1D)
177  {
178  s_D[k1][k2][k3] = D(k1,k2,k3,e);
179  }
180  }
181  }
182  MFEM_SYNC_THREAD;
183  MFEM_FOREACH_THREAD(i1,x,D1D)
184  {
185  MFEM_FOREACH_THREAD(i2,y,D1D)
186  {
187  MFEM_FOREACH_THREAD(i3,z,D1D)
188  {
189  for (int j1 = 0; j1 < D1D; ++j1)
190  {
191  for (int j2 = 0; j2 < D1D; ++j2)
192  {
193  for (int j3 = 0; j3 < D1D; ++j3)
194  {
195  double val = 0.0;
196  for (int k1 = 0; k1 < Q1D; ++k1)
197  {
198  for (int k2 = 0; k2 < Q1D; ++k2)
199  {
200  for (int k3 = 0; k3 < Q1D; ++k3)
201  {
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]
205  * s_D[k1][k2][k3];
206  }
207  }
208  }
209  if (add)
210  {
211  M(i1, i2, i3, j1, j2, j3, e) += val;
212  }
213  else
214  {
215  M(i1, i2, i3, j1, j2, j3, e) = val;
216  }
217  }
218  }
219  }
220  }
221  }
222  }
223  });
224 }
225 
226 void MassIntegrator::AssembleEA(const FiniteElementSpace &fes,
227  Vector &ea_data,
228  const bool add)
229 {
230  AssemblePA(fes);
231  const int ne = fes.GetMesh()->GetNE();
232  const Array<double> &B = maps->B;
233  if (dim == 1)
234  {
235  switch ((dofs1D << 4 ) | quad1D)
236  {
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,
246  dofs1D,quad1D);
247  }
248  }
249  else if (dim == 2)
250  {
251  switch ((dofs1D << 4 ) | quad1D)
252  {
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,
262  dofs1D,quad1D);
263  }
264  }
265  else if (dim == 3)
266  {
267  switch ((dofs1D << 4 ) | quad1D)
268  {
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,
277  dofs1D,quad1D);
278  }
279  }
280  MFEM_ABORT("Unknown kernel.");
281 }
282 
283 }
mfem::MassIntegrator::AssemblePA
virtual void AssemblePA(const FiniteElementSpace &fes)
Method defining partial assembly.
Definition: bilininteg_mass_pa.cpp:26
mfem::Mesh::GetNE
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:846
mfem::add
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:291
mfem::MassIntegrator::maps
const DofToQuad * maps
Not owned.
Definition: bilininteg.hpp:2066
mfem::Reshape
DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims...dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
Definition: dtensor.hpp:136
mfem::MAX_Q1D
const int MAX_Q1D
Definition: forall.hpp:28
mfem::FiniteElementSpace::GetMesh
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:410
mfem::MassIntegrator::AssembleEA
virtual void AssembleEA(const FiniteElementSpace &fes, Vector &emat, const bool add)
Method defining element assembly.
Definition: bilininteg_mass_ea.cpp:226
mfem::FiniteElementSpace
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
mfem::DofToQuad::B
Array< double > B
Basis functions evaluated at quadrature points.
Definition: fe.hpp:188
mfem::MAX_D1D
const int MAX_D1D
Definition: forall.hpp:27
mfem::Vector
Vector data type.
Definition: vector.hpp:60
Generated on Thu Jul 29 2021 15:58:37 for MFEM by   doxygen 1.8.5