MFEM  v4.2.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
bilininteg_dgtrace_ea.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, 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,
22  Vector &eadata_int,
23  Vector &eadata_ext,
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);
36  if (add)
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,
56  Vector &eadata_bdr,
57  const bool add)
58 {
59  auto D = Reshape(padata.Read(), 2, 2, NF);
60  auto A_bdr = Reshape(eadata_bdr.ReadWrite(), NF);
61  MFEM_FORALL(f, NF,
62  {
63  if (add)
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,
78  Vector &eadata_int,
79  Vector &eadata_ext,
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;
96  MFEM_FOREACH_THREAD(i1,x,D1D)
97  {
98  MFEM_FOREACH_THREAD(j1,y,D1D)
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  }
111  if (add)
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,
134  Vector &eadata_bdr,
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;
150  MFEM_FOREACH_THREAD(i1,x,D1D)
151  {
152  MFEM_FOREACH_THREAD(j1,y,D1D)
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  }
159  if (add)
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,
176  Vector &eadata_int,
177  Vector &eadata_ext,
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  {
209  MFEM_FOREACH_THREAD(k1,x,Q1D)
210  {
211  MFEM_FOREACH_THREAD(k2,y,Q1D)
212  {
213  s_D[k1][k2][i][j] = D(k1,k2,i,j,f);
214  }
215  }
216  }
217  }
218  MFEM_SYNC_THREAD;
219  MFEM_FOREACH_THREAD(i1,x,D1D)
220  {
221  MFEM_FOREACH_THREAD(i2,y,D1D)
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  }
249  if (add)
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,
274  Vector &eadata_bdr,
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  {
305  MFEM_FOREACH_THREAD(k1,x,Q1D)
306  {
307  MFEM_FOREACH_THREAD(k2,y,Q1D)
308  {
309  s_D[k1][k2][i][j] = D(k1,k2,i,j,f);
310  }
311  }
312  }
313  }
314  MFEM_SYNC_THREAD;
315  MFEM_FOREACH_THREAD(i1,x,D1D)
316  {
317  MFEM_FOREACH_THREAD(i2,y,D1D)
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  }
333  if (add)
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 
348 void DGTraceIntegrator::AssembleEAInteriorFaces(const FiniteElementSpace& fes,
349  Vector &ea_data_int,
350  Vector &ea_data_ext,
351  const bool add)
352 {
353  SetupPA(fes, FaceType::Interior);
354  nf = fes.GetNFbyType(FaceType::Interior);
355  if (nf==0) { return; }
356  const Array<double> &B = maps->B;
357  if (dim == 1)
358  {
359  return EADGTraceAssemble1DInt(nf,B,pa_data,ea_data_int,ea_data_ext,add);
360  }
361  else if (dim == 2)
362  {
363  switch ((dofs1D << 4 ) | quad1D)
364  {
365  case 0x22:
366  return EADGTraceAssemble2DInt<2,2>(nf,B,pa_data,ea_data_int,
367  ea_data_ext,add);
368  case 0x33:
369  return EADGTraceAssemble2DInt<3,3>(nf,B,pa_data,ea_data_int,
370  ea_data_ext,add);
371  case 0x44:
372  return EADGTraceAssemble2DInt<4,4>(nf,B,pa_data,ea_data_int,
373  ea_data_ext,add);
374  case 0x55:
375  return EADGTraceAssemble2DInt<5,5>(nf,B,pa_data,ea_data_int,
376  ea_data_ext,add);
377  case 0x66:
378  return EADGTraceAssemble2DInt<6,6>(nf,B,pa_data,ea_data_int,
379  ea_data_ext,add);
380  case 0x77:
381  return EADGTraceAssemble2DInt<7,7>(nf,B,pa_data,ea_data_int,
382  ea_data_ext,add);
383  case 0x88:
384  return EADGTraceAssemble2DInt<8,8>(nf,B,pa_data,ea_data_int,
385  ea_data_ext,add);
386  case 0x99:
387  return EADGTraceAssemble2DInt<9,9>(nf,B,pa_data,ea_data_int,
388  ea_data_ext,add);
389  default:
390  return EADGTraceAssemble2DInt(nf,B,pa_data,ea_data_int,
391  ea_data_ext,add,dofs1D,quad1D);
392  }
393  }
394  else if (dim == 3)
395  {
396  switch ((dofs1D << 4 ) | quad1D)
397  {
398  case 0x23:
399  return EADGTraceAssemble3DInt<2,3>(nf,B,pa_data,ea_data_int,
400  ea_data_ext,add);
401  case 0x34:
402  return EADGTraceAssemble3DInt<3,4>(nf,B,pa_data,ea_data_int,
403  ea_data_ext,add);
404  case 0x45:
405  return EADGTraceAssemble3DInt<4,5>(nf,B,pa_data,ea_data_int,
406  ea_data_ext,add);
407  case 0x56:
408  return EADGTraceAssemble3DInt<5,6>(nf,B,pa_data,ea_data_int,
409  ea_data_ext,add);
410  case 0x67:
411  return EADGTraceAssemble3DInt<6,7>(nf,B,pa_data,ea_data_int,
412  ea_data_ext,add);
413  case 0x78:
414  return EADGTraceAssemble3DInt<7,8>(nf,B,pa_data,ea_data_int,
415  ea_data_ext,add);
416  case 0x89:
417  return EADGTraceAssemble3DInt<8,9>(nf,B,pa_data,ea_data_int,
418  ea_data_ext,add);
419  default:
420  return EADGTraceAssemble3DInt(nf,B,pa_data,ea_data_int,
421  ea_data_ext,add,dofs1D,quad1D);
422  }
423  }
424  MFEM_ABORT("Unknown kernel.");
425 }
426 
427 void DGTraceIntegrator::AssembleEABoundaryFaces(const FiniteElementSpace& fes,
428  Vector &ea_data_bdr,
429  const bool add)
430 {
431  SetupPA(fes, FaceType::Boundary);
432  nf = fes.GetNFbyType(FaceType::Boundary);
433  if (nf==0) { return; }
434  const Array<double> &B = maps->B;
435  if (dim == 1)
436  {
437  return EADGTraceAssemble1DBdr(nf,B,pa_data,ea_data_bdr,add);
438  }
439  else if (dim == 2)
440  {
441  switch ((dofs1D << 4 ) | quad1D)
442  {
443  case 0x22: return EADGTraceAssemble2DBdr<2,2>(nf,B,pa_data,ea_data_bdr,add);
444  case 0x33: return EADGTraceAssemble2DBdr<3,3>(nf,B,pa_data,ea_data_bdr,add);
445  case 0x44: return EADGTraceAssemble2DBdr<4,4>(nf,B,pa_data,ea_data_bdr,add);
446  case 0x55: return EADGTraceAssemble2DBdr<5,5>(nf,B,pa_data,ea_data_bdr,add);
447  case 0x66: return EADGTraceAssemble2DBdr<6,6>(nf,B,pa_data,ea_data_bdr,add);
448  case 0x77: return EADGTraceAssemble2DBdr<7,7>(nf,B,pa_data,ea_data_bdr,add);
449  case 0x88: return EADGTraceAssemble2DBdr<8,8>(nf,B,pa_data,ea_data_bdr,add);
450  case 0x99: return EADGTraceAssemble2DBdr<9,9>(nf,B,pa_data,ea_data_bdr,add);
451  default:
452  return EADGTraceAssemble2DBdr(nf,B,pa_data,ea_data_bdr,add,dofs1D,quad1D);
453  }
454  }
455  else if (dim == 3)
456  {
457  switch ((dofs1D << 4 ) | quad1D)
458  {
459  case 0x23: return EADGTraceAssemble3DBdr<2,3>(nf,B,pa_data,ea_data_bdr,add);
460  case 0x34: return EADGTraceAssemble3DBdr<3,4>(nf,B,pa_data,ea_data_bdr,add);
461  case 0x45: return EADGTraceAssemble3DBdr<4,5>(nf,B,pa_data,ea_data_bdr,add);
462  case 0x56: return EADGTraceAssemble3DBdr<5,6>(nf,B,pa_data,ea_data_bdr,add);
463  case 0x67: return EADGTraceAssemble3DBdr<6,7>(nf,B,pa_data,ea_data_bdr,add);
464  case 0x78: return EADGTraceAssemble3DBdr<7,8>(nf,B,pa_data,ea_data_bdr,add);
465  case 0x89: return EADGTraceAssemble3DBdr<8,9>(nf,B,pa_data,ea_data_bdr,add);
466  default:
467  return EADGTraceAssemble3DBdr(nf,B,pa_data,ea_data_bdr,add,dofs1D,quad1D);
468  }
469  }
470  MFEM_ABORT("Unknown kernel.");
471 }
472 
473 }
mfem::FiniteElementSpace::GetNFbyType
int GetNFbyType(FaceType type) const
Returns the number of faces according to the requested type.
Definition: fespace.hpp:444
mfem::add
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:261
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:134
mfem::MAX_Q1D
const int MAX_Q1D
Definition: forall.hpp:27
mfem::FiniteElementSpace
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
mfem::DGTraceIntegrator::AssembleEAInteriorFaces
virtual void AssembleEAInteriorFaces(const FiniteElementSpace &fes, Vector &ea_data_int, Vector &ea_data_ext, const bool add)
Definition: bilininteg_dgtrace_ea.cpp:348
mfem::DGTraceIntegrator::AssembleEABoundaryFaces
virtual void AssembleEABoundaryFaces(const FiniteElementSpace &fes, Vector &ea_data_bdr, const bool add)
Definition: bilininteg_dgtrace_ea.cpp:427
mfem::DofToQuad::B
Array< double > B
Basis functions evaluated at quadrature points.
Definition: fe.hpp:180
mfem::MAX_D1D
const int MAX_D1D
Definition: forall.hpp:26
mfem::DGTraceIntegrator::maps
const DofToQuad * maps
Not owned.
Definition: bilininteg.hpp:2670
mfem::Vector
Vector data type.
Definition: vector.hpp:51
f
double f(const Vector &p)
Definition: mesh-explorer.cpp:76
Generated on Fri Oct 30 2020 15:47:59 for MFEM by   doxygen 1.8.5