MFEM  v4.5.2 Finite element discretization library
bilininteg_dgtrace_ea.cpp
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 #include "../general/forall.hpp"
13 #include "bilininteg.hpp"
14 #include "gridfunc.hpp"
15
16 namespace mfem
17 {
18
19 static void EADGTraceAssemble1DInt(const int NF,
20  const Array<double> &basis,
21  const Vector &padata,
24  const bool add)
25 {
26  auto D = Reshape(padata.Read(), 2, 2, NF);
27  auto A_int = Reshape(eadata_int.ReadWrite(), 2, NF);
28  auto A_ext = Reshape(eadata_ext.ReadWrite(), 2, NF);
29  MFEM_FORALL(f, NF,
30  {
31  double val_int0, val_int1, val_ext01, val_ext10;
32  val_int0 = D(0, 0, f);
33  val_ext10 = D(1, 0, f);
34  val_ext01 = D(0, 1, f);
35  val_int1 = D(1, 1, f);
37  {
38  A_int(0, f) += val_int0;
39  A_int(1, f) += val_int1;
40  A_ext(0, f) += val_ext01;
41  A_ext(1, f) += val_ext10;
42  }
43  else
44  {
45  A_int(0, f) = val_int0;
46  A_int(1, f) = val_int1;
47  A_ext(0, f) = val_ext01;
48  A_ext(1, f) = val_ext10;
49  }
50  });
51 }
52
53 static void EADGTraceAssemble1DBdr(const int NF,
54  const Array<double> &basis,
55  const Vector &padata,
57  const bool add)
58 {
59  auto D = Reshape(padata.Read(), 2, 2, NF);
61  MFEM_FORALL(f, NF,
62  {
64  {
65  A_bdr(f) += D(0, 0, f);
66  }
67  else
68  {
69  A_bdr(f) = D(0, 0, f);
70  }
71  });
72 }
73
74 template<int T_D1D = 0, int T_Q1D = 0>
75 static void EADGTraceAssemble2DInt(const int NF,
76  const Array<double> &basis,
77  const Vector &padata,
80  const bool add,
81  const int d1d = 0,
82  const int q1d = 0)
83 {
84  const int D1D = T_D1D ? T_D1D : d1d;
85  const int Q1D = T_Q1D ? T_Q1D : q1d;
86  MFEM_VERIFY(D1D <= MAX_D1D, "");
87  MFEM_VERIFY(Q1D <= MAX_Q1D, "");
88  auto B = Reshape(basis.Read(), Q1D, D1D);
89  auto D = Reshape(padata.Read(), Q1D, 2, 2, NF);
90  auto A_int = Reshape(eadata_int.ReadWrite(), D1D, D1D, 2, NF);
91  auto A_ext = Reshape(eadata_ext.ReadWrite(), D1D, D1D, 2, NF);
92  MFEM_FORALL_3D(f, NF, D1D, D1D, 1,
93  {
94  const int D1D = T_D1D ? T_D1D : d1d;
95  const int Q1D = T_Q1D ? T_Q1D : q1d;
97  {
99  {
100  double val_int0 = 0.0;
101  double val_int1 = 0.0;
102  double val_ext01 = 0.0;
103  double val_ext10 = 0.0;
104  for (int k1 = 0; k1 < Q1D; ++k1)
105  {
106  val_int0 += B(k1,i1) * B(k1,j1) * D(k1, 0, 0, f);
107  val_ext01 += B(k1,i1) * B(k1,j1) * D(k1, 0, 1, f);
108  val_ext10 += B(k1,i1) * B(k1,j1) * D(k1, 1, 0, f);
109  val_int1 += B(k1,i1) * B(k1,j1) * D(k1, 1, 1, f);
110  }
112  {
113  A_int(i1, j1, 0, f) += val_int0;
114  A_int(i1, j1, 1, f) += val_int1;
115  A_ext(i1, j1, 0, f) += val_ext01;
116  A_ext(i1, j1, 1, f) += val_ext10;
117  }
118  else
119  {
120  A_int(i1, j1, 0, f) = val_int0;
121  A_int(i1, j1, 1, f) = val_int1;
122  A_ext(i1, j1, 0, f) = val_ext01;
123  A_ext(i1, j1, 1, f) = val_ext10;
124  }
125  }
126  }
127  });
128 }
129
130 template<int T_D1D = 0, int T_Q1D = 0>
131 static void EADGTraceAssemble2DBdr(const int NF,
132  const Array<double> &basis,
133  const Vector &padata,
135  const bool add,
136  const int d1d = 0,
137  const int q1d = 0)
138 {
139  const int D1D = T_D1D ? T_D1D : d1d;
140  const int Q1D = T_Q1D ? T_Q1D : q1d;
141  MFEM_VERIFY(D1D <= MAX_D1D, "");
142  MFEM_VERIFY(Q1D <= MAX_Q1D, "");
143  auto B = Reshape(basis.Read(), Q1D, D1D);
144  auto D = Reshape(padata.Read(), Q1D, 2, 2, NF);
145  auto A_bdr = Reshape(eadata_bdr.ReadWrite(), D1D, D1D, NF);
146  MFEM_FORALL_3D(f, NF, D1D, D1D, 1,
147  {
148  const int D1D = T_D1D ? T_D1D : d1d;
149  const int Q1D = T_Q1D ? T_Q1D : q1d;
151  {
153  {
154  double val_bdr = 0.0;
155  for (int k1 = 0; k1 < Q1D; ++k1)
156  {
157  val_bdr += B(k1,i1) * B(k1,j1) * D(k1, 0, 0, f);
158  }
160  {
161  A_bdr(i1, j1, f) += val_bdr;
162  }
163  else
164  {
165  A_bdr(i1, j1, f) = val_bdr;
166  }
167  }
168  }
169  });
170 }
171
172 template<int T_D1D = 0, int T_Q1D = 0>
173 static void EADGTraceAssemble3DInt(const int NF,
174  const Array<double> &basis,
175  const Vector &padata,
178  const bool add,
179  const int d1d = 0,
180  const int q1d = 0)
181 {
182  const int D1D = T_D1D ? T_D1D : d1d;
183  const int Q1D = T_Q1D ? T_Q1D : q1d;
184  MFEM_VERIFY(D1D <= MAX_D1D, "");
185  MFEM_VERIFY(Q1D <= MAX_Q1D, "");
186  auto B = Reshape(basis.Read(), Q1D, D1D);
187  auto D = Reshape(padata.Read(), Q1D, Q1D, 2, 2, NF);
188  auto A_int = Reshape(eadata_int.ReadWrite(), D1D, D1D, D1D, D1D, 2, NF);
189  auto A_ext = Reshape(eadata_ext.ReadWrite(), D1D, D1D, D1D, D1D, 2, NF);
190  MFEM_FORALL_3D(f, NF, D1D, D1D, 1,
191  {
192  const int D1D = T_D1D ? T_D1D : d1d;
193  const int Q1D = T_Q1D ? T_Q1D : q1d;
194  constexpr int MD1 = T_D1D ? T_D1D : MAX_D1D;
195  constexpr int MQ1 = T_Q1D ? T_Q1D : MAX_Q1D;
196  double r_B[MQ1][MD1];
197  for (int d = 0; d < D1D; d++)
198  {
199  for (int q = 0; q < Q1D; q++)
200  {
201  r_B[q][d] = B(q,d);
202  }
203  }
204  MFEM_SHARED double s_D[MQ1][MQ1][2][2];
205  for (int i=0; i < 2; i++)
206  {
207  for (int j=0; j < 2; j++)
208  {
210  {
212  {
213  s_D[k1][k2][i][j] = D(k1,k2,i,j,f);
214  }
215  }
216  }
217  }
220  {
222  {
223  for (int j1 = 0; j1 < D1D; ++j1)
224  {
225  for (int j2 = 0; j2 < D1D; ++j2)
226  {
227  double val_int0 = 0.0;
228  double val_int1 = 0.0;
229  double val_ext01 = 0.0;
230  double val_ext10 = 0.0;
231  for (int k1 = 0; k1 < Q1D; ++k1)
232  {
233  for (int k2 = 0; k2 < Q1D; ++k2)
234  {
235  val_int0 += r_B[k1][i1] * r_B[k1][j1]
236  * r_B[k2][i2] * r_B[k2][j2]
237  * s_D[k1][k2][0][0];
238  val_int1 += r_B[k1][i1] * r_B[k1][j1]
239  * r_B[k2][i2] * r_B[k2][j2]
240  * s_D[k1][k2][1][1];
241  val_ext01+= r_B[k1][i1] * r_B[k1][j1]
242  * r_B[k2][i2] * r_B[k2][j2]
243  * s_D[k1][k2][0][1];
244  val_ext10+= r_B[k1][i1] * r_B[k1][j1]
245  * r_B[k2][i2] * r_B[k2][j2]
246  * s_D[k1][k2][1][0];
247  }
248  }
250  {
251  A_int(i1, i2, j1, j2, 0, f) += val_int0;
252  A_int(i1, i2, j1, j2, 1, f) += val_int1;
253  A_ext(i1, i2, j1, j2, 0, f) += val_ext01;
254  A_ext(i1, i2, j1, j2, 1, f) += val_ext10;
255  }
256  else
257  {
258  A_int(i1, i2, j1, j2, 0, f) = val_int0;
259  A_int(i1, i2, j1, j2, 1, f) = val_int1;
260  A_ext(i1, i2, j1, j2, 0, f) = val_ext01;
261  A_ext(i1, i2, j1, j2, 1, f) = val_ext10;
262  }
263  }
264  }
265  }
266  }
267  });
268 }
269
270 template<int T_D1D = 0, int T_Q1D = 0>
271 static void EADGTraceAssemble3DBdr(const int NF,
272  const Array<double> &basis,
273  const Vector &padata,
275  const bool add,
276  const int d1d = 0,
277  const int q1d = 0)
278 {
279  const int D1D = T_D1D ? T_D1D : d1d;
280  const int Q1D = T_Q1D ? T_Q1D : q1d;
281  MFEM_VERIFY(D1D <= MAX_D1D, "");
282  MFEM_VERIFY(Q1D <= MAX_Q1D, "");
283  auto B = Reshape(basis.Read(), Q1D, D1D);
284  auto D = Reshape(padata.Read(), Q1D, Q1D, 2, 2, NF);
285  auto A_bdr = Reshape(eadata_bdr.ReadWrite(), D1D, D1D, D1D, D1D, NF);
286  MFEM_FORALL_3D(f, NF, D1D, D1D, 1,
287  {
288  const int D1D = T_D1D ? T_D1D : d1d;
289  const int Q1D = T_Q1D ? T_Q1D : q1d;
290  constexpr int MD1 = T_D1D ? T_D1D : MAX_D1D;
291  constexpr int MQ1 = T_Q1D ? T_Q1D : MAX_Q1D;
292  double r_B[MQ1][MD1];
293  for (int d = 0; d < D1D; d++)
294  {
295  for (int q = 0; q < Q1D; q++)
296  {
297  r_B[q][d] = B(q,d);
298  }
299  }
300  MFEM_SHARED double s_D[MQ1][MQ1][2][2];
301  for (int i=0; i < 2; i++)
302  {
303  for (int j=0; j < 2; j++)
304  {
306  {
308  {
309  s_D[k1][k2][i][j] = D(k1,k2,i,j,f);
310  }
311  }
312  }
313  }
316  {
318  {
319  for (int j1 = 0; j1 < D1D; ++j1)
320  {
321  for (int j2 = 0; j2 < D1D; ++j2)
322  {
323  double val_bdr = 0.0;
324  for (int k1 = 0; k1 < Q1D; ++k1)
325  {
326  for (int k2 = 0; k2 < Q1D; ++k2)
327  {
328  val_bdr += r_B[k1][i1] * r_B[k1][j1]
329  * r_B[k2][i2] * r_B[k2][j2]
330  * s_D[k1][k2][0][0];
331  }
332  }
334  {
335  A_bdr(i1, i2, j1, j2, f) += val_bdr;
336  }
337  else
338  {
339  A_bdr(i1, i2, j1, j2, f) = val_bdr;
340  }
341  }
342  }
343  }
344  }
345  });
346 }
347
349  Vector &ea_data_int,
350  Vector &ea_data_ext,
351  const bool add)
352 {
353  SetupPA(fes, FaceType::Interior);
355  if (nf==0) { return; }
356  const Array<double> &B = maps->B;
357  if (dim == 1)
358  {
360  }
361  else if (dim == 2)
362  {
363  switch ((dofs1D << 4 ) | quad1D)
364  {
365  case 0x22:
368  case 0x33:
371  case 0x44:
374  case 0x55:
377  case 0x66:
380  case 0x77:
383  case 0x88:
386  case 0x99:
389  default:
392  }
393  }
394  else if (dim == 3)
395  {
396  switch ((dofs1D << 4 ) | quad1D)
397  {
398  case 0x23:
401  case 0x34:
404  case 0x45:
407  case 0x56:
410  case 0x67:
413  case 0x78:
416  case 0x89:
419  default:
422  }
423  }
424  MFEM_ABORT("Unknown kernel.");
425 }
426
428  Vector &ea_data_bdr,
429  const bool add)
430 {
431  SetupPA(fes, FaceType::Boundary);
433  if (nf==0) { return; }
434  const Array<double> &B = maps->B;
435  if (dim == 1)
436  {
438  }
439  else if (dim == 2)
440  {
441  switch ((dofs1D << 4 ) | quad1D)
442  {
451  default:
453  }
454  }
455  else if (dim == 3)
456  {
457  switch ((dofs1D << 4 ) | quad1D)
458  {
466  default:
468  }
469  }
470  MFEM_ABORT("Unknown kernel.");
471 }
472
473 }
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:314
int GetNFbyType(FaceType type) const
Returns the number of faces according to the requested type.
Definition: fespace.hpp:631
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
const int MAX_Q1D
Definition: forall.hpp:29
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
virtual void AssembleEAInteriorFaces(const FiniteElementSpace &fes, Vector &ea_data_int, Vector &ea_data_ext, const bool add)
virtual void AssembleEABoundaryFaces(const FiniteElementSpace &fes, Vector &ea_data_bdr, const bool add)
Array< double > B
Basis functions evaluated at quadrature points.
Definition: fe_base.hpp:183
const int MAX_D1D
Definition: forall.hpp:28
const DofToQuad * maps
Not owned.
Vector data type.
Definition: vector.hpp:60
MFEM_HOST_DEVICE DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims... dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
Definition: dtensor.hpp:131