MFEM  v4.2.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
bilininteg_vecmass.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 #include "libceed/mass.hpp"
16 
17 using namespace std;
18 
19 namespace mfem
20 {
21 
22 // PA Mass Integrator
23 
24 // PA Mass Assemble kernel
25 void VectorMassIntegrator::AssemblePA(const FiniteElementSpace &fes)
26 {
27  // Assuming the same element type
28  Mesh *mesh = fes.GetMesh();
29  if (mesh->GetNE() == 0) { return; }
30  const FiniteElement &el = *fes.GetFE(0);
31  ElementTransformation *T = mesh->GetElementTransformation(0);
32  const IntegrationRule *ir
33  = IntRule ? IntRule : &MassIntegrator::GetRule(el, el, *T);
34  if (DeviceCanUseCeed())
35  {
36  delete ceedDataPtr;
37  ceedDataPtr = new CeedData;
38  InitCeedCoeff(Q, *mesh, *ir, ceedDataPtr);
39  return CeedPAMassAssemble(fes, *ir, *ceedDataPtr);
40  }
41  dim = mesh->Dimension();
42  ne = fes.GetMesh()->GetNE();
43  nq = ir->GetNPoints();
44  geom = mesh->GetGeometricFactors(*ir, GeometricFactors::COORDINATES |
45  GeometricFactors::JACOBIANS);
46  maps = &el.GetDofToQuad(*ir, DofToQuad::TENSOR);
47  dofs1D = maps->ndof;
48  quad1D = maps->nqpt;
49  pa_data.SetSize(ne*nq, Device::GetDeviceMemoryType());
50  double coeff = 1.0;
51  if (Q)
52  {
53  ConstantCoefficient *cQ = dynamic_cast<ConstantCoefficient*>(Q);
54  MFEM_VERIFY(cQ != NULL, "Only ConstantCoefficient is supported.");
55  coeff = cQ->constant;
56  }
57  if (!(dim == 2 || dim == 3))
58  {
59  MFEM_ABORT("Dimension not supported.");
60  }
61  if (dim == 2)
62  {
63  const double constant = coeff;
64  const int NE = ne;
65  const int NQ = nq;
66  auto w = ir->GetWeights().Read();
67  auto J = Reshape(geom->J.Read(), NQ,2,2,NE);
68  auto v = Reshape(pa_data.Write(), NQ, NE);
69  MFEM_FORALL(e, NE,
70  {
71  for (int q = 0; q < NQ; ++q)
72  {
73  const double J11 = J(q,0,0,e);
74  const double J12 = J(q,1,0,e);
75  const double J21 = J(q,0,1,e);
76  const double J22 = J(q,1,1,e);
77  const double detJ = (J11*J22)-(J21*J12);
78  v(q,e) = w[q] * constant * detJ;
79  }
80  });
81  }
82  if (dim == 3)
83  {
84  const double constant = coeff;
85  const int NE = ne;
86  const int NQ = nq;
87  auto W = ir->GetWeights().Read();
88  auto J = Reshape(geom->J.Read(), NQ,3,3,NE);
89  auto v = Reshape(pa_data.Write(), NQ,NE);
90  MFEM_FORALL(e, NE,
91  {
92  for (int q = 0; q < NQ; ++q)
93  {
94  const double J11 = J(q,0,0,e), J12 = J(q,0,1,e), J13 = J(q,0,2,e);
95  const double J21 = J(q,1,0,e), J22 = J(q,1,1,e), J23 = J(q,1,2,e);
96  const double J31 = J(q,2,0,e), J32 = J(q,2,1,e), J33 = J(q,2,2,e);
97  const double detJ = J11 * (J22 * J33 - J32 * J23) -
98  /* */ J21 * (J12 * J33 - J32 * J13) +
99  /* */ J31 * (J12 * J23 - J22 * J13);
100  v(q,e) = W[q] * constant * detJ;
101  }
102  });
103  }
104 }
105 
106 template<const int T_D1D = 0,
107  const int T_Q1D = 0>
108 static void PAVectorMassApply2D(const int NE,
109  const Array<double> &_B,
110  const Array<double> &_Bt,
111  const Vector &_op,
112  const Vector &_x,
113  Vector &_y,
114  const int d1d = 0,
115  const int q1d = 0)
116 {
117  const int D1D = T_D1D ? T_D1D : d1d;
118  const int Q1D = T_Q1D ? T_Q1D : q1d;
119  constexpr int VDIM = 2;
120  MFEM_VERIFY(D1D <= MAX_D1D, "");
121  MFEM_VERIFY(Q1D <= MAX_Q1D, "");
122  auto B = Reshape(_B.Read(), Q1D, D1D);
123  auto Bt = Reshape(_Bt.Read(), D1D, Q1D);
124  auto op = Reshape(_op.Read(), Q1D, Q1D, NE);
125  auto x = Reshape(_x.Read(), D1D, D1D, VDIM, NE);
126  auto y = Reshape(_y.ReadWrite(), D1D, D1D, VDIM, NE);
127  MFEM_FORALL(e, NE,
128  {
129  const int D1D = T_D1D ? T_D1D : d1d; // nvcc workaround
130  const int Q1D = T_Q1D ? T_Q1D : q1d;
131  // the following variables are evaluated at compile time
132  constexpr int max_D1D = T_D1D ? T_D1D : MAX_D1D;
133  constexpr int max_Q1D = T_Q1D ? T_Q1D : MAX_Q1D;
134  double sol_xy[max_Q1D][max_Q1D];
135  for (int c = 0; c < VDIM; ++c)
136  {
137  for (int qy = 0; qy < Q1D; ++qy)
138  {
139  for (int qx = 0; qx < Q1D; ++qx)
140  {
141  sol_xy[qy][qx] = 0.0;
142  }
143  }
144  for (int dy = 0; dy < D1D; ++dy)
145  {
146  double sol_x[max_Q1D];
147  for (int qy = 0; qy < Q1D; ++qy)
148  {
149  sol_x[qy] = 0.0;
150  }
151  for (int dx = 0; dx < D1D; ++dx)
152  {
153  const double s = x(dx,dy,c,e);
154  for (int qx = 0; qx < Q1D; ++qx)
155  {
156  sol_x[qx] += B(qx,dx)* s;
157  }
158  }
159  for (int qy = 0; qy < Q1D; ++qy)
160  {
161  const double d2q = B(qy,dy);
162  for (int qx = 0; qx < Q1D; ++qx)
163  {
164  sol_xy[qy][qx] += d2q * sol_x[qx];
165  }
166  }
167  }
168  for (int qy = 0; qy < Q1D; ++qy)
169  {
170  for (int qx = 0; qx < Q1D; ++qx)
171  {
172  sol_xy[qy][qx] *= op(qx,qy,e);
173  }
174  }
175  for (int qy = 0; qy < Q1D; ++qy)
176  {
177  double sol_x[max_D1D];
178  for (int dx = 0; dx < D1D; ++dx)
179  {
180  sol_x[dx] = 0.0;
181  }
182  for (int qx = 0; qx < Q1D; ++qx)
183  {
184  const double s = sol_xy[qy][qx];
185  for (int dx = 0; dx < D1D; ++dx)
186  {
187  sol_x[dx] += Bt(dx,qx) * s;
188  }
189  }
190  for (int dy = 0; dy < D1D; ++dy)
191  {
192  const double q2d = Bt(dy,qy);
193  for (int dx = 0; dx < D1D; ++dx)
194  {
195  y(dx,dy,c,e) += q2d * sol_x[dx];
196  }
197  }
198  }
199  }
200  });
201 }
202 
203 template<const int T_D1D = 0,
204  const int T_Q1D = 0>
205 static void PAVectorMassApply3D(const int NE,
206  const Array<double> &_B,
207  const Array<double> &_Bt,
208  const Vector &_op,
209  const Vector &_x,
210  Vector &_y,
211  const int d1d = 0,
212  const int q1d = 0)
213 {
214  const int D1D = T_D1D ? T_D1D : d1d;
215  const int Q1D = T_Q1D ? T_Q1D : q1d;
216  constexpr int VDIM = 3;
217  MFEM_VERIFY(D1D <= MAX_D1D, "");
218  MFEM_VERIFY(Q1D <= MAX_Q1D, "");
219  auto B = Reshape(_B.Read(), Q1D, D1D);
220  auto Bt = Reshape(_Bt.Read(), D1D, Q1D);
221  auto op = Reshape(_op.Read(), Q1D, Q1D, Q1D, NE);
222  auto x = Reshape(_x.Read(), D1D, D1D, D1D, VDIM, NE);
223  auto y = Reshape(_y.ReadWrite(), D1D, D1D, D1D, VDIM, NE);
224  MFEM_FORALL(e, NE,
225  {
226  const int D1D = T_D1D ? T_D1D : d1d;
227  const int Q1D = T_Q1D ? T_Q1D : q1d;
228  constexpr int max_D1D = T_D1D ? T_D1D : MAX_D1D;
229  constexpr int max_Q1D = T_Q1D ? T_Q1D : MAX_Q1D;
230  double sol_xyz[max_Q1D][max_Q1D][max_Q1D];
231  for (int c = 0; c < VDIM; ++ c)
232  {
233  for (int qz = 0; qz < Q1D; ++qz)
234  {
235  for (int qy = 0; qy < Q1D; ++qy)
236  {
237  for (int qx = 0; qx < Q1D; ++qx)
238  {
239  sol_xyz[qz][qy][qx] = 0.0;
240  }
241  }
242  }
243  for (int dz = 0; dz < D1D; ++dz)
244  {
245  double sol_xy[max_Q1D][max_Q1D];
246  for (int qy = 0; qy < Q1D; ++qy)
247  {
248  for (int qx = 0; qx < Q1D; ++qx)
249  {
250  sol_xy[qy][qx] = 0.0;
251  }
252  }
253  for (int dy = 0; dy < D1D; ++dy)
254  {
255  double sol_x[max_Q1D];
256  for (int qx = 0; qx < Q1D; ++qx)
257  {
258  sol_x[qx] = 0;
259  }
260  for (int dx = 0; dx < D1D; ++dx)
261  {
262  const double s = x(dx,dy,dz,c,e);
263  for (int qx = 0; qx < Q1D; ++qx)
264  {
265  sol_x[qx] += B(qx,dx) * s;
266  }
267  }
268  for (int qy = 0; qy < Q1D; ++qy)
269  {
270  const double wy = B(qy,dy);
271  for (int qx = 0; qx < Q1D; ++qx)
272  {
273  sol_xy[qy][qx] += wy * sol_x[qx];
274  }
275  }
276  }
277  for (int qz = 0; qz < Q1D; ++qz)
278  {
279  const double wz = B(qz,dz);
280  for (int qy = 0; qy < Q1D; ++qy)
281  {
282  for (int qx = 0; qx < Q1D; ++qx)
283  {
284  sol_xyz[qz][qy][qx] += wz * sol_xy[qy][qx];
285  }
286  }
287  }
288  }
289  for (int qz = 0; qz < Q1D; ++qz)
290  {
291  for (int qy = 0; qy < Q1D; ++qy)
292  {
293  for (int qx = 0; qx < Q1D; ++qx)
294  {
295  sol_xyz[qz][qy][qx] *= op(qx,qy,qz,e);
296  }
297  }
298  }
299  for (int qz = 0; qz < Q1D; ++qz)
300  {
301  double sol_xy[max_D1D][max_D1D];
302  for (int dy = 0; dy < D1D; ++dy)
303  {
304  for (int dx = 0; dx < D1D; ++dx)
305  {
306  sol_xy[dy][dx] = 0;
307  }
308  }
309  for (int qy = 0; qy < Q1D; ++qy)
310  {
311  double sol_x[max_D1D];
312  for (int dx = 0; dx < D1D; ++dx)
313  {
314  sol_x[dx] = 0;
315  }
316  for (int qx = 0; qx < Q1D; ++qx)
317  {
318  const double s = sol_xyz[qz][qy][qx];
319  for (int dx = 0; dx < D1D; ++dx)
320  {
321  sol_x[dx] += Bt(dx,qx) * s;
322  }
323  }
324  for (int dy = 0; dy < D1D; ++dy)
325  {
326  const double wy = Bt(dy,qy);
327  for (int dx = 0; dx < D1D; ++dx)
328  {
329  sol_xy[dy][dx] += wy * sol_x[dx];
330  }
331  }
332  }
333  for (int dz = 0; dz < D1D; ++dz)
334  {
335  const double wz = Bt(dz,qz);
336  for (int dy = 0; dy < D1D; ++dy)
337  {
338  for (int dx = 0; dx < D1D; ++dx)
339  {
340  y(dx,dy,dz,c,e) += wz * sol_xy[dy][dx];
341  }
342  }
343  }
344  }
345  }
346  });
347 }
348 
349 static void PAVectorMassApply(const int dim,
350  const int D1D,
351  const int Q1D,
352  const int NE,
353  const Array<double> &B,
354  const Array<double> &Bt,
355  const Vector &op,
356  const Vector &x,
357  Vector &y)
358 {
359  if (dim == 2)
360  {
361  return PAVectorMassApply2D(NE, B, Bt, op, x, y, D1D, Q1D);
362  }
363  if (dim == 3)
364  {
365  return PAVectorMassApply3D(NE, B, Bt, op, x, y, D1D, Q1D);
366  }
367  MFEM_ABORT("Unknown kernel.");
368 }
369 
370 void VectorMassIntegrator::AddMultPA(const Vector &x, Vector &y) const
371 {
372  if (DeviceCanUseCeed())
373  {
374  CeedAddMult(ceedDataPtr, x, y);
375  }
376  else
377  {
378  PAVectorMassApply(dim, dofs1D, quad1D, ne, maps->B, maps->Bt, pa_data, x, y);
379  }
380 }
381 
382 template<const int T_D1D = 0, const int T_Q1D = 0>
383 static void PAVectorMassAssembleDiagonal2D(const int NE,
384  const Array<double> &_B,
385  const Array<double> &_Bt,
386  const Vector &_op,
387  Vector &_diag,
388  const int d1d = 0,
389  const int q1d = 0)
390 {
391  const int D1D = T_D1D ? T_D1D : d1d;
392  const int Q1D = T_Q1D ? T_Q1D : q1d;
393  constexpr int VDIM = 2;
394  MFEM_VERIFY(D1D <= MAX_D1D, "");
395  MFEM_VERIFY(Q1D <= MAX_Q1D, "");
396  auto B = Reshape(_B.Read(), Q1D, D1D);
397  auto op = Reshape(_op.Read(), Q1D, Q1D, NE);
398  auto y = Reshape(_diag.ReadWrite(), D1D, D1D, VDIM, NE);
399  MFEM_FORALL(e, NE,
400  {
401  const int D1D = T_D1D ? T_D1D : d1d;
402  const int Q1D = T_Q1D ? T_Q1D : q1d;
403  constexpr int max_D1D = T_D1D ? T_D1D : MAX_D1D;
404  constexpr int max_Q1D = T_Q1D ? T_Q1D : MAX_Q1D;
405 
406  double temp[max_Q1D][max_D1D];
407  for (int qx = 0; qx < Q1D; ++qx)
408  {
409  for (int dy = 0; dy < D1D; ++dy)
410  {
411  temp[qx][dy] = 0.0;
412  for (int qy = 0; qy < Q1D; ++qy)
413  {
414  temp[qx][dy] += B(qy, dy) * B(qy, dy) * op(qx, qy, e);
415  }
416  }
417  }
418  for (int dy = 0; dy < D1D; ++dy)
419  {
420  for (int dx = 0; dx < D1D; ++dx)
421  {
422  double temp1 = 0.0;
423  for (int qx = 0; qx < Q1D; ++qx)
424  {
425  temp1 += B(qx, dx) * B(qx, dx) * temp[qx][dy];
426  }
427  y(dx, dy, 0, e) = temp1;
428  y(dx, dy, 1, e) = temp1;
429  }
430  }
431  });
432 }
433 
434 template<const int T_D1D = 0, const int T_Q1D = 0>
435 static void PAVectorMassAssembleDiagonal3D(const int NE,
436  const Array<double> &_B,
437  const Array<double> &_Bt,
438  const Vector &_op,
439  Vector &_diag,
440  const int d1d = 0,
441  const int q1d = 0)
442 {
443  const int D1D = T_D1D ? T_D1D : d1d;
444  const int Q1D = T_Q1D ? T_Q1D : q1d;
445  constexpr int VDIM = 3;
446  MFEM_VERIFY(D1D <= MAX_D1D, "");
447  MFEM_VERIFY(Q1D <= MAX_Q1D, "");
448  auto B = Reshape(_B.Read(), Q1D, D1D);
449  auto op = Reshape(_op.Read(), Q1D, Q1D, Q1D, NE);
450  auto y = Reshape(_diag.ReadWrite(), D1D, D1D, D1D, VDIM, NE);
451  MFEM_FORALL(e, NE,
452  {
453  const int D1D = T_D1D ? T_D1D : d1d; // nvcc workaround
454  const int Q1D = T_Q1D ? T_Q1D : q1d;
455  // the following variables are evaluated at compile time
456  constexpr int max_D1D = T_D1D ? T_D1D : MAX_D1D;
457  constexpr int max_Q1D = T_Q1D ? T_Q1D : MAX_Q1D;
458 
459  double temp[max_Q1D][max_Q1D][max_D1D];
460  for (int qx = 0; qx < Q1D; ++qx)
461  {
462  for (int qy = 0; qy < Q1D; ++qy)
463  {
464  for (int dz = 0; dz < D1D; ++dz)
465  {
466  temp[qx][qy][dz] = 0.0;
467  for (int qz = 0; qz < Q1D; ++qz)
468  {
469  temp[qx][qy][dz] += B(qz, dz) * B(qz, dz) * op(qx, qy, qz, e);
470  }
471  }
472  }
473  }
474  double temp2[max_Q1D][max_D1D][max_D1D];
475  for (int qx = 0; qx < Q1D; ++qx)
476  {
477  for (int dz = 0; dz < D1D; ++dz)
478  {
479  for (int dy = 0; dy < D1D; ++dy)
480  {
481  temp2[qx][dy][dz] = 0.0;
482  for (int qy = 0; qy < Q1D; ++qy)
483  {
484  temp2[qx][dy][dz] += B(qy, dy) * B(qy, dy) * temp[qx][qy][dz];
485  }
486  }
487  }
488  }
489  for (int dz = 0; dz < D1D; ++dz)
490  {
491  for (int dy = 0; dy < D1D; ++dy)
492  {
493  for (int dx = 0; dx < D1D; ++dx)
494  {
495  double temp3 = 0.0;
496  for (int qx = 0; qx < Q1D; ++qx)
497  {
498  temp3 += B(qx, dx) * B(qx, dx)
499  * temp2[qx][dy][dz];
500  }
501  y(dx, dy, dz, 0, e) = temp3;
502  y(dx, dy, dz, 1, e) = temp3;
503  y(dx, dy, dz, 2, e) = temp3;
504  }
505  }
506  }
507  });
508 }
509 
510 static void PAVectorMassAssembleDiagonal(const int dim,
511  const int D1D,
512  const int Q1D,
513  const int NE,
514  const Array<double> &B,
515  const Array<double> &Bt,
516  const Vector &op,
517  Vector &y)
518 {
519  if (dim == 2)
520  {
521  return PAVectorMassAssembleDiagonal2D(NE, B, Bt, op, y, D1D, Q1D);
522  }
523  else if (dim == 3)
524  {
525  return PAVectorMassAssembleDiagonal3D(NE, B, Bt, op, y, D1D, Q1D);
526  }
527  MFEM_ABORT("Dimension not implemented.");
528 }
529 
530 void VectorMassIntegrator::AssembleDiagonalPA(Vector &diag)
531 {
532  if (DeviceCanUseCeed())
533  {
534  CeedAssembleDiagonal(ceedDataPtr, diag);
535  }
536  else
537  {
538  PAVectorMassAssembleDiagonal(dim,
539  dofs1D,
540  quad1D,
541  ne,
542  maps->B,
543  maps->Bt,
544  pa_data,
545  diag);
546  }
547 }
548 
549 } // namespace mfem
mfem::IntegrationRule::GetNPoints
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:245
mfem::FiniteElement
Abstract class for all finite elements.
Definition: fe.hpp:235
mfem::IntegrationRule
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
mfem::ConstantCoefficient
A coefficient that is constant across space and time.
Definition: coefficient.hpp:78
geom
const Geometry::Type geom
Definition: ex1.cpp:40
mfem::DofToQuad::ndof
int ndof
Number of degrees of freedom = number of basis functions. When mode is TENSOR, this is the 1D number...
Definition: fe.hpp:165
mfem::Mesh::GetNE
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:737
mfem::IntegrationRule::GetWeights
const Array< double > & GetWeights() const
Return the quadrature weights in a contiguous array.
Definition: intrules.cpp:85
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::Mesh::GetGeometricFactors
const GeometricFactors * GetGeometricFactors(const IntegrationRule &ir, const int flags)
Return the mesh geometric factors corresponding to the given integration rule.
Definition: mesh.cpp:768
mfem::Vector::ReadWrite
double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:388
mfem::MAX_Q1D
const int MAX_Q1D
Definition: forall.hpp:27
mfem::FiniteElementSpace::GetMesh
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:314
mfem::Array::Read
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:290
mfem::Mesh::Dimension
int Dimension() const
Definition: mesh.hpp:788
mfem::FiniteElementSpace
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
mfem::FiniteElement::GetDofToQuad
virtual const DofToQuad & GetDofToQuad(const IntegrationRule &ir, DofToQuad::Mode mode) const
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
Definition: fe.cpp:376
mfem::Vector::Read
const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:372
mfem::Mesh::GetElementTransformation
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:336
dim
int dim
Definition: ex24.cpp:53
mfem::MAX_D1D
const int MAX_D1D
Definition: forall.hpp:26
mfem::FiniteElementSpace::GetFE
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:1798
mfem::Vector
Vector data type.
Definition: vector.hpp:51
Generated on Fri Oct 30 2020 15:47:59 for MFEM by   doxygen 1.8.5