MFEM  v4.5.2
Finite element discretization library
lor_h1.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 "lor_h1.hpp"
13 #include "lor_util.hpp"
14 #include "../../linalg/dtensor.hpp"
15 #include "../../general/forall.hpp"
16 
17 namespace mfem
18 {
19 
20 template <int ORDER>
22 {
23  const int nel_ho = fes_ho.GetNE();
24 
25  static constexpr int nv = 4;
26  static constexpr int dim = 2;
27  static constexpr int ddm2 = (dim*(dim+1))/2;
28  static constexpr int nd1d = ORDER + 1;
29  static constexpr int ndof_per_el = nd1d*nd1d;
30  static constexpr int nnz_per_row = 9;
31  static constexpr int sz_local_mat = nv*nv;
32 
33  const bool const_mq = c1.Size() == 1;
34  const auto MQ = const_mq
35  ? Reshape(c1.Read(), 1, 1, 1)
36  : Reshape(c1.Read(), nd1d, nd1d, nel_ho);
37  const bool const_dq = c2.Size() == 1;
38  const auto DQ = const_dq
39  ? Reshape(c2.Read(), 1, 1, 1)
40  : Reshape(c2.Read(), nd1d, nd1d, nel_ho);
41 
42  sparse_ij.SetSize(nnz_per_row*ndof_per_el*nel_ho);
43  auto V = Reshape(sparse_ij.Write(), nnz_per_row, nd1d, nd1d, nel_ho);
44 
45  auto X = X_vert.Read();
46 
47  MFEM_FORALL_2D(iel_ho, nel_ho, ORDER, ORDER, 1,
48  {
49  // Assemble a sparse matrix over the macro-element by looping over each
50  // subelement.
51  // V(j,ix,iy) stores the jth nonzero in the row of the sparse matrix
52  // corresponding to local DOF (ix, iy).
53  MFEM_FOREACH_THREAD(iy,y,nd1d)
54  {
55  MFEM_FOREACH_THREAD(ix,x,nd1d)
56  {
57  for (int j=0; j<nnz_per_row; ++j)
58  {
59  V(j,ix,iy,iel_ho) = 0.0;
60  }
61  }
62  }
63  MFEM_SYNC_THREAD;
64 
65  // Compute geometric factors at quadrature points
66  MFEM_FOREACH_THREAD(ky,y,ORDER)
67  {
68  MFEM_FOREACH_THREAD(kx,x,ORDER)
69  {
70  double Q_[(ddm2 + 1)*nv];
71  double local_mat_[sz_local_mat];
72  DeviceTensor<3> Q(Q_, ddm2 + 1, 2, 2);
73  DeviceTensor<2> local_mat(local_mat_, nv, nv);
74 
75  for (int i=0; i<sz_local_mat; ++i) { local_mat[i] = 0.0; }
76 
77  double vx[4], vy[4];
78  LORVertexCoordinates2D<ORDER>(X, iel_ho, kx, ky, vx, vy);
79 
80  for (int iqy=0; iqy<2; ++iqy)
81  {
82  for (int iqx=0; iqx<2; ++iqx)
83  {
84  const double x = iqx;
85  const double y = iqy;
86  const double w = 1.0/4.0;
87 
88  double J_[2*2];
89  DeviceTensor<2> J(J_, 2, 2);
90 
91  Jacobian2D(x, y, vx, vy, J);
92 
93  const double detJ = Det2D(J);
94  const double w_detJ = w/detJ;
95 
96  Q(0,iqy,iqx) = w_detJ * (J(0,1)*J(0,1) + J(1,1)*J(1,1)); // 1,1
97  Q(1,iqy,iqx) = -w_detJ * (J(0,1)*J(0,0) + J(1,1)*J(1,0)); // 1,2
98  Q(2,iqy,iqx) = w_detJ * (J(0,0)*J(0,0) + J(1,0)*J(1,0)); // 2,2
99  Q(3,iqy,iqx) = w*detJ;
100  }
101  }
102  for (int iqx=0; iqx<2; ++iqx)
103  {
104  for (int iqy=0; iqy<2; ++iqy)
105  {
106  const double mq = const_mq ? MQ(0,0,0) : MQ(kx+iqx, ky+iqy, iel_ho);
107  const double dq = const_dq ? DQ(0,0,0) : DQ(kx+iqx, ky+iqy, iel_ho);
108  for (int jy=0; jy<2; ++jy)
109  {
110  const double bjy = (jy == iqy) ? 1.0 : 0.0;
111  const double gjy = (jy == 0) ? -1.0 : 1.0;
112  for (int jx=0; jx<2; ++jx)
113  {
114  const double bjx = (jx == iqx) ? 1.0 : 0.0;
115  const double gjx = (jx == 0) ? -1.0 : 1.0;
116 
117  const double djx = gjx*bjy;
118  const double djy = bjx*gjy;
119 
120  int jj_loc = jx + 2*jy;
121 
122  for (int iy=0; iy<2; ++iy)
123  {
124  const double biy = (iy == iqy) ? 1.0 : 0.0;
125  const double giy = (iy == 0) ? -1.0 : 1.0;
126  for (int ix=0; ix<2; ++ix)
127  {
128  const double bix = (ix == iqx) ? 1.0 : 0.0;
129  const double gix = (ix == 0) ? -1.0 : 1.0;
130 
131  const double dix = gix*biy;
132  const double diy = bix*giy;
133 
134  int ii_loc = ix + 2*iy;
135 
136  // Only store the lower-triangular part of
137  // the matrix (by symmetry).
138  if (jj_loc > ii_loc) { continue; }
139 
140  double val = 0.0;
141  val += dix*djx*Q(0,iqy,iqx);
142  val += (dix*djy + diy*djx)*Q(1,iqy,iqx);
143  val += diy*djy*Q(2,iqy,iqx);
144  val *= dq;
145 
146  val += mq*bix*biy*bjx*bjy*Q(3,iqy,iqx);
147 
148  local_mat(ii_loc, jj_loc) += val;
149  }
150  }
151  }
152  }
153  }
154  }
155  // Assemble the local matrix into the macro-element sparse matrix
156  // in a format similar to coordinate format. The (I,J) arrays
157  // are implicit (not stored explicitly).
158  for (int ii_loc=0; ii_loc<nv; ++ii_loc)
159  {
160  const int ix = ii_loc%2;
161  const int iy = ii_loc/2;
162  for (int jj_loc=0; jj_loc<nv; ++jj_loc)
163  {
164  const int jx = jj_loc%2;
165  const int jy = jj_loc/2;
166  const int jj_off = (jx-ix+1) + 3*(jy-iy+1);
167 
168  // Symmetry
169  if (jj_loc <= ii_loc)
170  {
171  AtomicAdd(V(jj_off, ix+kx, iy+ky, iel_ho), local_mat(ii_loc, jj_loc));
172  }
173  else
174  {
175  AtomicAdd(V(jj_off, ix+kx, iy+ky, iel_ho), local_mat(jj_loc, ii_loc));
176  }
177  }
178  }
179  }
180  }
181  });
182 
183  sparse_mapping.SetSize(nnz_per_row*ndof_per_el);
184  sparse_mapping = -1;
185  auto map = Reshape(sparse_mapping.HostReadWrite(), nnz_per_row, ndof_per_el);
186  for (int iy=0; iy<nd1d; ++iy)
187  {
188  const int jy_begin = (iy > 0) ? iy - 1 : 0;
189  const int jy_end = (iy < ORDER) ? iy + 1 : ORDER;
190  for (int ix=0; ix<nd1d; ++ix)
191  {
192  const int jx_begin = (ix > 0) ? ix - 1 : 0;
193  const int jx_end = (ix < ORDER) ? ix + 1 : ORDER;
194  const int ii_el = ix + nd1d*iy;
195  for (int jy=jy_begin; jy<=jy_end; ++jy)
196  {
197  for (int jx=jx_begin; jx<=jx_end; ++jx)
198  {
199  const int jj_off = (jx-ix+1) + 3*(jy-iy+1);
200  const int jj_el = jx + nd1d*jy;
201  map(jj_off, ii_el) = jj_el;
202  }
203  }
204  }
205  }
206 }
207 
208 template <int ORDER>
210 {
211  const int nel_ho = fes_ho.GetNE();
212  static constexpr int nv = 8;
213  static constexpr int dim = 3;
214  static constexpr int ddm2 = (dim*(dim+1))/2;
215  static constexpr int nd1d = ORDER + 1;
216  static constexpr int ndof_per_el = nd1d*nd1d*nd1d;
217  static constexpr int nnz_per_row = 27;
218  static constexpr int sz_grad_A = 3*3*2*2*2*2;
219  static constexpr int sz_grad_B = sz_grad_A*2;
220  static constexpr int sz_mass_A = 2*2*2*2;
221  static constexpr int sz_mass_B = sz_mass_A*2;
222  static constexpr int sz_local_mat = nv*nv;
223 
224  const bool const_mq = c1.Size() == 1;
225  const auto MQ = const_mq
226  ? Reshape(c1.Read(), 1, 1, 1, 1)
227  : Reshape(c1.Read(), nd1d, nd1d, nd1d, nel_ho);
228  const bool const_dq = c2.Size() == 1;
229  const auto DQ = const_dq
230  ? Reshape(c2.Read(), 1, 1, 1, 1)
231  : Reshape(c2.Read(), nd1d, nd1d, nd1d, nel_ho);
232 
233  sparse_ij.SetSize(nel_ho*ndof_per_el*nnz_per_row);
234  auto V = Reshape(sparse_ij.Write(), nnz_per_row, nd1d, nd1d, nd1d, nel_ho);
235 
236  auto X = X_vert.Read();
237 
238  // Last thread dimension is lowered to avoid "too many resources" error
239  MFEM_FORALL_3D(iel_ho, nel_ho, ORDER, ORDER, (ORDER>6)?4:ORDER,
240  {
241  // Assemble a sparse matrix over the macro-element by looping over each
242  // subelement.
243  // V(j,i) stores the jth nonzero in the ith row of the sparse matrix.
244  MFEM_FOREACH_THREAD(iz,z,nd1d)
245  {
246  MFEM_FOREACH_THREAD(iy,y,nd1d)
247  {
248  MFEM_FOREACH_THREAD(ix,x,nd1d)
249  {
250  MFEM_UNROLL(nnz_per_row)
251  for (int j=0; j<nnz_per_row; ++j)
252  {
253  V(j,ix,iy,iz,iel_ho) = 0.0;
254  }
255  }
256  }
257  }
258  MFEM_SYNC_THREAD;
259 
260  // Compute geometric factors at quadrature points
261  MFEM_FOREACH_THREAD(kz,z,ORDER)
262  {
263  MFEM_FOREACH_THREAD(ky,y,ORDER)
264  {
265  MFEM_FOREACH_THREAD(kx,x,ORDER)
266  {
267  double Q_[(ddm2 + 1)*nv];
268  double grad_A_[sz_grad_A];
269  double grad_B_[sz_grad_B];
270  double mass_A_[sz_mass_A];
271  double mass_B_[sz_mass_B];
272  double local_mat_[sz_local_mat];
273 
274  DeviceTensor<4> Q(Q_, ddm2 + 1, 2, 2, 2);
275  DeviceTensor<2> local_mat(local_mat_, nv, nv);
276  DeviceTensor<6> grad_A(grad_A_, 3, 3, 2, 2, 2, 2);
277  DeviceTensor<7> grad_B(grad_B_, 3, 3, 2, 2, 2, 2, 2);
278  DeviceTensor<4> mass_A(mass_A_, 2, 2, 2, 2);
279  DeviceTensor<5> mass_B(mass_B_, 2, 2, 2, 2, 2);
280 
281  // local_mat is the local (dense) stiffness matrix
282  for (int i=0; i<sz_local_mat; ++i) { local_mat[i] = 0.0; }
283 
284  // Intermediate quantities
285  // (see e.g. Mora and Demkowicz for notation).
286  for (int i=0; i<sz_grad_A; ++i) { grad_A[i] = 0.0; }
287  for (int i=0; i<sz_grad_B; ++i) { grad_B[i] = 0.0; }
288 
289  for (int i=0; i<sz_mass_A; ++i) { mass_A[i] = 0.0; }
290  for (int i=0; i<sz_mass_B; ++i) { mass_B[i] = 0.0; }
291 
292  double vx[8], vy[8], vz[8];
293  LORVertexCoordinates3D<ORDER>(X, iel_ho, kx, ky, kz, vx, vy, vz);
294 
295  //MFEM_UNROLL(2)
296  for (int iqz=0; iqz<2; ++iqz)
297  {
298  //MFEM_UNROLL(2)
299  for (int iqy=0; iqy<2; ++iqy)
300  {
301  //MFEM_UNROLL(2)
302  for (int iqx=0; iqx<2; ++iqx)
303  {
304  const double x = iqx;
305  const double y = iqy;
306  const double z = iqz;
307  const double w = 1.0/8.0;
308 
309  double J_[3*3];
310  DeviceTensor<2> J(J_, 3, 3);
311 
312  Jacobian3D(x, y, z, vx, vy, vz, J);
313 
314  const double detJ = Det3D(J);
315  const double w_detJ = w/detJ;
316 
317  // adj(J)
318  double A_[3*3];
319  DeviceTensor<2> A(A_, 3, 3);
320  Adjugate3D(J, A);
321 
322  Q(0,iqz,iqy,iqx) = w_detJ*(A(0,0)*A(0,0)+A(0,1)*A(0,1)+A(0,2)*A(0,2)); // 1,1
323  Q(1,iqz,iqy,iqx) = w_detJ*(A(0,0)*A(1,0)+A(0,1)*A(1,1)+A(0,2)*A(1,2)); // 2,1
324  Q(2,iqz,iqy,iqx) = w_detJ*(A(0,0)*A(2,0)+A(0,1)*A(2,1)+A(0,2)*A(2,2)); // 3,1
325  Q(3,iqz,iqy,iqx) = w_detJ*(A(1,0)*A(1,0)+A(1,1)*A(1,1)+A(1,2)*A(1,2)); // 2,2
326  Q(4,iqz,iqy,iqx) = w_detJ*(A(1,0)*A(2,0)+A(1,1)*A(2,1)+A(1,2)*A(2,2)); // 3,2
327  Q(5,iqz,iqy,iqx) = w_detJ*(A(2,0)*A(2,0)+A(2,1)*A(2,1)+A(2,2)*A(2,2)); // 3,3
328  Q(6,iqz,iqy,iqx) = w*detJ;
329  }
330  }
331  }
332 
333  //MFEM_UNROLL(2)
334  for (int iqx=0; iqx<2; ++iqx)
335  {
336  //MFEM_UNROLL(2)
337  for (int jz=0; jz<2; ++jz)
338  {
339  // Note loop starts at iz=jz here, taking advantage of
340  // symmetries.
341  //MFEM_UNROLL(2)
342  for (int iz=jz; iz<2; ++iz)
343  {
344  //MFEM_UNROLL(2)
345  for (int iqy=0; iqy<2; ++iqy)
346  {
347  //MFEM_UNROLL(2)
348  for (int iqz=0; iqz<2; ++iqz)
349  {
350  const double mq = const_mq ? MQ(0,0,0,0) : MQ(kx+iqx, ky+iqy, kz+iqz, iel_ho);
351  const double dq = const_dq ? DQ(0,0,0,0) : DQ(kx+iqx, ky+iqy, kz+iqz, iel_ho);
352 
353  const double biz = (iz == iqz) ? 1.0 : 0.0;
354  const double giz = (iz == 0) ? -1.0 : 1.0;
355 
356  const double bjz = (jz == iqz) ? 1.0 : 0.0;
357  const double gjz = (jz == 0) ? -1.0 : 1.0;
358 
359  const double J11 = Q(0,iqz,iqy,iqx);
360  const double J21 = Q(1,iqz,iqy,iqx);
361  const double J31 = Q(2,iqz,iqy,iqx);
362  const double J12 = J21;
363  const double J22 = Q(3,iqz,iqy,iqx);
364  const double J32 = Q(4,iqz,iqy,iqx);
365  const double J13 = J31;
366  const double J23 = J32;
367  const double J33 = Q(5,iqz,iqy,iqx);
368 
369  grad_A(0,0,iqy,iz,jz,iqx) += dq*J11*biz*bjz;
370  grad_A(1,0,iqy,iz,jz,iqx) += dq*J21*biz*bjz;
371  grad_A(2,0,iqy,iz,jz,iqx) += dq*J31*giz*bjz;
372  grad_A(0,1,iqy,iz,jz,iqx) += dq*J12*biz*bjz;
373  grad_A(1,1,iqy,iz,jz,iqx) += dq*J22*biz*bjz;
374  grad_A(2,1,iqy,iz,jz,iqx) += dq*J32*giz*bjz;
375  grad_A(0,2,iqy,iz,jz,iqx) += dq*J13*biz*gjz;
376  grad_A(1,2,iqy,iz,jz,iqx) += dq*J23*biz*gjz;
377  grad_A(2,2,iqy,iz,jz,iqx) += dq*J33*giz*gjz;
378 
379  double wdetJ = Q(6,iqz,iqy,iqx);
380  mass_A(iqy,iz,jz,iqx) += mq*wdetJ*biz*bjz;
381  }
382  //MFEM_UNROLL(2)
383  for (int jy=0; jy<2; ++jy)
384  {
385  //MFEM_UNROLL(2)
386  for (int iy=0; iy<2; ++iy)
387  {
388  const double biy = (iy == iqy) ? 1.0 : 0.0;
389  const double giy = (iy == 0) ? -1.0 : 1.0;
390 
391  const double bjy = (jy == iqy) ? 1.0 : 0.0;
392  const double gjy = (jy == 0) ? -1.0 : 1.0;
393 
394  grad_B(0,0,iy,jy,iz,jz,iqx) += biy*bjy*grad_A(0,0,iqy,iz,jz,iqx);
395  grad_B(1,0,iy,jy,iz,jz,iqx) += giy*bjy*grad_A(1,0,iqy,iz,jz,iqx);
396  grad_B(2,0,iy,jy,iz,jz,iqx) += biy*bjy*grad_A(2,0,iqy,iz,jz,iqx);
397  grad_B(0,1,iy,jy,iz,jz,iqx) += biy*gjy*grad_A(0,1,iqy,iz,jz,iqx);
398  grad_B(1,1,iy,jy,iz,jz,iqx) += giy*gjy*grad_A(1,1,iqy,iz,jz,iqx);
399  grad_B(2,1,iy,jy,iz,jz,iqx) += biy*gjy*grad_A(2,1,iqy,iz,jz,iqx);
400  grad_B(0,2,iy,jy,iz,jz,iqx) += biy*bjy*grad_A(0,2,iqy,iz,jz,iqx);
401  grad_B(1,2,iy,jy,iz,jz,iqx) += giy*bjy*grad_A(1,2,iqy,iz,jz,iqx);
402  grad_B(2,2,iy,jy,iz,jz,iqx) += biy*bjy*grad_A(2,2,iqy,iz,jz,iqx);
403 
404  mass_B(iy,jy,iz,jz,iqx) += biy*bjy*mass_A(iqy,iz,jz,iqx);
405  }
406  }
407  }
408  //MFEM_UNROLL(2)
409  for (int jy=0; jy<2; ++jy)
410  {
411  //MFEM_UNROLL(2)
412  for (int jx=0; jx<2; ++jx)
413  {
414  //MFEM_UNROLL(2)
415  for (int iy=0; iy<2; ++iy)
416  {
417  //MFEM_UNROLL(2)
418  for (int ix=0; ix<2; ++ix)
419  {
420  const double bix = (ix == iqx) ? 1.0 : 0.0;
421  const double gix = (ix == 0) ? -1.0 : 1.0;
422 
423  const double bjx = (jx == iqx) ? 1.0 : 0.0;
424  const double gjx = (jx == 0) ? -1.0 : 1.0;
425 
426  int ii_loc = ix + 2*iy + 4*iz;
427  int jj_loc = jx + 2*jy + 4*jz;
428 
429  // Only store the lower-triangular part of
430  // the matrix (by symmetry).
431  if (jj_loc > ii_loc) { continue; }
432 
433  double val = 0.0;
434  val += gix*gjx*grad_B(0,0,iy,jy,iz,jz,iqx);
435  val += bix*gjx*grad_B(1,0,iy,jy,iz,jz,iqx);
436  val += bix*gjx*grad_B(2,0,iy,jy,iz,jz,iqx);
437  val += gix*bjx*grad_B(0,1,iy,jy,iz,jz,iqx);
438  val += bix*bjx*grad_B(1,1,iy,jy,iz,jz,iqx);
439  val += bix*bjx*grad_B(2,1,iy,jy,iz,jz,iqx);
440  val += gix*bjx*grad_B(0,2,iy,jy,iz,jz,iqx);
441  val += bix*bjx*grad_B(2,2,iy,jy,iz,jz,iqx);
442  val += bix*bjx*grad_B(1,2,iy,jy,iz,jz,iqx);
443 
444  val += bix*bjx*mass_B(iy,jy,iz,jz,iqx);
445 
446  local_mat(ii_loc, jj_loc) += val;
447  }
448  }
449  }
450  }
451  }
452  }
453  }
454  // Assemble the local matrix into the macro-element sparse matrix
455  // in a format similar to coordinate format. The (I,J) arrays
456  // are implicit (not stored explicitly).
457  //MFEM_UNROLL(8)
458  for (int ii_loc=0; ii_loc<nv; ++ii_loc)
459  {
460  const int ix = ii_loc%2;
461  const int iy = (ii_loc/2)%2;
462  const int iz = ii_loc/2/2;
463 
464  for (int jj_loc=0; jj_loc<nv; ++jj_loc)
465  {
466  const int jx = jj_loc%2;
467  const int jy = (jj_loc/2)%2;
468  const int jz = jj_loc/2/2;
469  const int jj_off = (jx-ix+1) + 3*(jy-iy+1) + 9*(jz-iz+1);
470 
471  if (jj_loc <= ii_loc)
472  {
473  AtomicAdd(V(jj_off, ix+kx, iy+ky, iz+kz, iel_ho), local_mat(ii_loc, jj_loc));
474  }
475  else
476  {
477  AtomicAdd(V(jj_off, ix+kx, iy+ky, iz+kz, iel_ho), local_mat(jj_loc, ii_loc));
478  }
479  }
480  }
481  }
482  }
483  }
484  });
485 
486  sparse_mapping.SetSize(nnz_per_row*ndof_per_el);
487  sparse_mapping = -1;
488  auto map = Reshape(sparse_mapping.HostReadWrite(), nnz_per_row, ndof_per_el);
489  for (int iz=0; iz<nd1d; ++iz)
490  {
491  const int jz_begin = (iz > 0) ? iz - 1 : 0;
492  const int jz_end = (iz < ORDER) ? iz + 1 : ORDER;
493  for (int iy=0; iy<nd1d; ++iy)
494  {
495  const int jy_begin = (iy > 0) ? iy - 1 : 0;
496  const int jy_end = (iy < ORDER) ? iy + 1 : ORDER;
497  for (int ix=0; ix<nd1d; ++ix)
498  {
499  const int jx_begin = (ix > 0) ? ix - 1 : 0;
500  const int jx_end = (ix < ORDER) ? ix + 1 : ORDER;
501 
502  const int ii_el = ix + nd1d*(iy + nd1d*iz);
503 
504  for (int jz=jz_begin; jz<=jz_end; ++jz)
505  {
506  for (int jy=jy_begin; jy<=jy_end; ++jy)
507  {
508  for (int jx=jx_begin; jx<=jx_end; ++jx)
509  {
510  const int jj_off = (jx-ix+1) + 3*(jy-iy+1) + 9*(jz-iz+1);
511  const int jj_el = jx + nd1d*(jy + nd1d*jz);
512  map(jj_off, ii_el) = jj_el;
513  }
514  }
515  }
516  }
517  }
518  }
519 }
520 
521 // Explicit template instantiations
522 template void BatchedLOR_H1::Assemble2D<1>();
523 template void BatchedLOR_H1::Assemble2D<2>();
524 template void BatchedLOR_H1::Assemble2D<3>();
525 template void BatchedLOR_H1::Assemble2D<4>();
526 template void BatchedLOR_H1::Assemble2D<5>();
527 template void BatchedLOR_H1::Assemble2D<6>();
528 template void BatchedLOR_H1::Assemble2D<7>();
529 template void BatchedLOR_H1::Assemble2D<8>();
530 
531 template void BatchedLOR_H1::Assemble3D<1>();
532 template void BatchedLOR_H1::Assemble3D<2>();
533 template void BatchedLOR_H1::Assemble3D<3>();
534 template void BatchedLOR_H1::Assemble3D<4>();
535 template void BatchedLOR_H1::Assemble3D<5>();
536 template void BatchedLOR_H1::Assemble3D<6>();
537 template void BatchedLOR_H1::Assemble3D<7>();
538 template void BatchedLOR_H1::Assemble3D<8>();
539 
541  FiniteElementSpace &fes_ho_,
542  Vector &X_vert_,
543  Vector &sparse_ij_,
544  Array<int> &sparse_mapping_)
545  : BatchedLORKernel(fes_ho_, X_vert_, sparse_ij_, sparse_mapping_)
546 {
547  ProjectLORCoefficient<MassIntegrator>(a, c1);
548  ProjectLORCoefficient<DiffusionIntegrator>(a, c2);
549 }
550 
551 } // namespace mfem
T * HostReadWrite()
Shortcut for mfem::ReadWrite(a.GetMemory(), a.Size(), false).
Definition: array.hpp:327
Abstract base class for the batched LOR assembly kernels.
Array< int > & sparse_mapping
Local element sparsity pattern.
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:512
MFEM_HOST_DEVICE T AtomicAdd(T &add, const T val)
Definition: backends.hpp:84
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:448
CoefficientVector c2
Coefficient of second integrator.
virtual double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:456
Vector & sparse_ij
Local element sparsity matrix data.
MFEM_HOST_DEVICE double Det3D(DeviceMatrix &J)
Definition: lor_util.hpp:188
MFEM_HOST_DEVICE double Det2D(DeviceMatrix &J)
Definition: lor_util.hpp:183
MFEM_HOST_DEVICE void Jacobian3D(const double x, const double y, const double z, const double vx[8], const double vy[8], const double vz[8], DeviceMatrix &J)
Definition: lor_util.hpp:126
FiniteElementSpace & fes_ho
The associated high-order space.
A basic generic Tensor class, appropriate for use on the GPU.
Definition: dtensor.hpp:81
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:614
MFEM_HOST_DEVICE void Jacobian2D(const double x, const double y, const double vx[4], const double vy[4], DeviceMatrix &J)
Definition: lor_util.hpp:115
Vector & X_vert
Mesh coordinate vector.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:684
BatchedLOR_H1(BilinearForm &a, FiniteElementSpace &fes_ho_, Vector &X_vert_, Vector &sparse_ij_, Array< int > &sparse_mapping_)
Definition: lor_h1.cpp:540
double a
Definition: lissajous.cpp:41
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
int dim
Definition: ex24.cpp:53
Vector data type.
Definition: vector.hpp:60
MFEM_HOST_DEVICE void Adjugate3D(const DeviceMatrix &J, DeviceMatrix &A)
Definition: lor_util.hpp:170
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
CoefficientVector c1
Coefficient of first integrator.