MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
bilininteg_dgtrace_ea.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, 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
16namespace mfem
17{
18
19static void EADGTraceAssemble1DInt(const int NF,
20 const Array<real_t> &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(NF, [=] MFEM_HOST_DEVICE (int f)
30 {
31 real_t 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
53static void EADGTraceAssemble1DBdr(const int NF,
54 const Array<real_t> &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(NF, [=] MFEM_HOST_DEVICE (int f)
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
74template<int T_D1D = 0, int T_Q1D = 0>
75static void EADGTraceAssemble2DInt(const int NF,
76 const Array<real_t> &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 <= DeviceDofQuadLimits::Get().MAX_D1D, "");
87 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().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_2D(NF, D1D, D1D, [=] MFEM_HOST_DEVICE (int f)
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 real_t val_int0 = 0.0;
101 real_t val_int1 = 0.0;
102 real_t val_ext01 = 0.0;
103 real_t 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
130template<int T_D1D = 0, int T_Q1D = 0>
131static void EADGTraceAssemble2DBdr(const int NF,
132 const Array<real_t> &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 <= DeviceDofQuadLimits::Get().MAX_D1D, "");
142 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().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_2D(NF, D1D, D1D, [=] MFEM_HOST_DEVICE (int f)
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 real_t 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
172template<int T_D1D = 0, int T_Q1D = 0>
173static void EADGTraceAssemble3DInt(const int NF,
174 const Array<real_t> &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 <= DeviceDofQuadLimits::Get().MAX_D1D, "");
185 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().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_2D(NF, D1D, D1D, [=] MFEM_HOST_DEVICE (int f)
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 : DofQuadLimits::MAX_D1D;
195 constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
196 real_t 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 real_t 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 real_t val_int0 = 0.0;
228 real_t val_int1 = 0.0;
229 real_t val_ext01 = 0.0;
230 real_t 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
270template<int T_D1D = 0, int T_Q1D = 0>
271static void EADGTraceAssemble3DBdr(const int NF,
272 const Array<real_t> &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 <= DeviceDofQuadLimits::Get().MAX_D1D, "");
282 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().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_2D(NF, D1D, D1D, [=] MFEM_HOST_DEVICE (int f)
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 : DofQuadLimits::MAX_D1D;
291 constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
292 real_t 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 real_t 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 real_t 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
348void 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<real_t> &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
427void 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<real_t> &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::DGTraceIntegrator::AssembleEABoundaryFaces
void AssembleEABoundaryFaces(const FiniteElementSpace &fes, Vector &ea_data_bdr, const bool add) override
Definition bilininteg_dgtrace_ea.cpp:427
mfem::DGTraceIntegrator::AssembleEAInteriorFaces
void AssembleEAInteriorFaces(const FiniteElementSpace &fes, Vector &ea_data_int, Vector &ea_data_ext, const bool add) override
Definition bilininteg_dgtrace_ea.cpp:348
mfem::DGTraceIntegrator::maps
const DofToQuad * maps
Not owned.
Definition bilininteg.hpp:3297
mfem::DofToQuad::B
Array< real_t > B
Basis functions evaluated at quadrature points.
Definition fe_base.hpp:193
mfem::FiniteElementSpace
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:244
mfem::FiniteElementSpace::GetNFbyType
int GetNFbyType(FaceType type) const
Returns the number of faces according to the requested type.
Definition fespace.hpp:908
mfem::Vector
Vector data type.
Definition vector.hpp:82
mfem::add
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition vector.cpp:391
mfem::Reshape
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
mfem::forall_2D
void forall_2D(int N, int X, int Y, lambda &&body)
Definition forall.hpp:762
mfem::real_t
float real_t
Definition config.hpp:43
mfem::f
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
mfem::forall
void forall(int N, lambda &&body)
Definition forall.hpp:753
mfem::DeviceDofQuadLimits::Get
static const DeviceDofQuadLimits & Get()
Return a const reference to the DeviceDofQuadLimits singleton.
Definition forall.hpp:128
Generated on Wed Apr 9 2025 16:53:04 for MFEM by doxygen 1.11.0