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