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