MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
lininteg_domain_grad.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 "fem.hpp"
13 #include "../fem/kernels.hpp"
14 #include "../general/forall.hpp"
15 
16 namespace mfem
17 {
18 
19 template<int T_D1D = 0, int T_Q1D = 0> static
20 void DLFGradAssemble2D(const int vdim, const int ne, const int d, const int q,
21  const int *markers, const double *b, const double *g,
22  const double *jacobians,
23  const double *weights, const Vector &coeff, double *y)
24 {
25  const auto F = coeff.Read();
26  const auto M = Reshape(markers, ne);
27  const auto B = Reshape(b, q, d);
28  const auto G = Reshape(g, q, d);
29  const auto J = Reshape(jacobians, q, q, 2,2, ne);
30  const auto W = Reshape(weights, q, q);
31  const bool cst = coeff.Size() == vdim*2;
32  const auto C = cst ? Reshape(F,2,vdim,1,1,1) : Reshape(F,2,vdim,q,q,ne);
33  auto Y = Reshape(y, d,d, vdim, ne);
34 
35  MFEM_FORALL_2D(e, ne, q, q, 1,
36  {
37  if (M(e) == 0) { return; } // ignore
38 
39  constexpr int Q = T_Q1D ? T_Q1D : MAX_Q1D;
40  constexpr int D = T_D1D ? T_D1D : MAX_D1D;
41 
42  MFEM_SHARED double sBGt[2][Q*D];
43  MFEM_SHARED double sQQ[2][Q*Q];
44  MFEM_SHARED double sDQ[2][D*Q];
45 
46  const DeviceMatrix Bt(sBGt[0], q, d);
47  const DeviceMatrix Gt(sBGt[1], q, d);
48  kernels::internal::LoadBGt<D,Q>(d, q, B, G, sBGt);
49 
50  const DeviceMatrix QQ0(sQQ[0], q, q);
51  const DeviceMatrix QQ1(sQQ[1], q, q);
52 
53  const DeviceMatrix DQ0(sDQ[0], d, q);
54  const DeviceMatrix DQ1(sDQ[1], d, q);
55 
56  for (int c = 0; c < vdim; ++c)
57  {
58  const double cst_val0 = C(0,c,0,0,0);
59  const double cst_val1 = C(1,c,0,0,0);
60 
61  MFEM_FOREACH_THREAD(x,x,q)
62  {
63  MFEM_FOREACH_THREAD(y,y,q)
64  {
65  const double w = W(x,y);
66  const double J11 = J(x,y,0,0,e);
67  const double J21 = J(x,y,1,0,e);
68  const double J12 = J(x,y,0,1,e);
69  const double J22 = J(x,y,1,1,e);
70  const double u = cst ? cst_val0 : C(0,c,x,y,e);
71  const double v = cst ? cst_val1 : C(1,c,x,y,e);
72  // QQ = w * det(J) * J^{-1} . C = w * adj(J) . { u, v }
73  QQ0(y,x) = w * (J22*u - J12*v);
74  QQ1(y,x) = w * (J11*v - J21*u);
75  }
76  }
77  MFEM_SYNC_THREAD;
78  MFEM_FOREACH_THREAD(qx,x,q)
79  {
80  MFEM_FOREACH_THREAD(dy,y,d)
81  {
82  double u = 0.0, v = 0.0;
83  for (int qy = 0; qy < q; ++qy)
84  {
85  u += QQ0(qy,qx) * Bt(qy,dy);
86  v += QQ1(qy,qx) * Gt(qy,dy);
87  }
88  DQ0(dy,qx) = u;
89  DQ1(dy,qx) = v;
90  }
91  }
92  MFEM_SYNC_THREAD;
93  MFEM_FOREACH_THREAD(dx,x,d)
94  {
95  MFEM_FOREACH_THREAD(dy,y,d)
96  {
97  double u = 0.0, v = 0.0;
98  for (int qx = 0; qx < q; ++qx)
99  {
100  u += DQ0(dy,qx) * Gt(qx,dx);
101  v += DQ1(dy,qx) * Bt(qx,dx);
102  }
103  Y(dx,dy,c,e) += u + v;
104  }
105  }
106  MFEM_SYNC_THREAD;
107  }
108  });
109 }
110 
111 template<int T_D1D = 0, int T_Q1D = 0> static
112 void DLFGradAssemble3D(const int vdim, const int ne, const int d, const int q,
113  const int *markers, const double *b, const double *g,
114  const double *jacobians,
115  const double *weights, const Vector &coeff,
116  double *output)
117 {
118  const auto F = coeff.Read();
119  const auto M = Reshape(markers, ne);
120  const auto B = Reshape(b, q,d);
121  const auto G = Reshape(g, q,d);
122  const auto J = Reshape(jacobians, q,q,q, 3,3, ne);
123  const auto W = Reshape(weights, q,q,q);
124  const bool cst = coeff.Size() == vdim*3;
125  const auto C = cst ? Reshape(F,3,vdim,1,1,1,1) : Reshape(F,3,vdim,q,q,q,ne);
126 
127  auto Y = Reshape(output, d,d,d, vdim, ne);
128 
129  MFEM_FORALL_2D(e, ne, q, q, 1,
130  {
131  if (M(e) == 0) { return; } // ignore
132 
133  constexpr int Q = T_Q1D ? T_Q1D : MAX_Q1D;
134  constexpr int D = T_D1D ? T_D1D : MAX_D1D;
135 
136  double r_u[D];
137 
138  MFEM_SHARED double sBGt[2][Q*D];
139  MFEM_SHARED double sQQQ[Q*Q*Q];
140 
141  const DeviceMatrix Bt(sBGt[0], q,d), Gt(sBGt[1], q,d);
142  kernels::internal::LoadBGt<D,Q>(d,q,B,G,sBGt);
143 
144  const DeviceCube QQQ(sQQQ, q,q,q);
145  const DeviceCube QQD(sQQQ, q,q,d);
146  const DeviceCube QDD(sQQQ, q,d,d);
147 
148  for (int c = 0; c < vdim; ++c)
149  {
150  const double cst_val_0 = C(0,c,0,0,0,0);
151  const double cst_val_1 = C(1,c,0,0,0,0);
152  const double cst_val_2 = C(2,c,0,0,0,0);
153 
154  for (int k = 0; k < 3; ++k)
155  {
156  for (int z = 0; z < q; ++z)
157  {
158  MFEM_FOREACH_THREAD(y,y,q)
159  {
160  MFEM_FOREACH_THREAD(x,x,q)
161  {
162  const double J11 = J(x,y,z,0,0,e);
163  const double J21 = J(x,y,z,1,0,e);
164  const double J31 = J(x,y,z,2,0,e);
165  const double J12 = J(x,y,z,0,1,e);
166  const double J22 = J(x,y,z,1,1,e);
167  const double J32 = J(x,y,z,2,1,e);
168  const double J13 = J(x,y,z,0,2,e);
169  const double J23 = J(x,y,z,1,2,e);
170  const double J33 = J(x,y,z,2,2,e);
171 
172  const double u = cst ? cst_val_0 : C(0,c,x,y,z,e);
173  const double v = cst ? cst_val_1 : C(1,c,x,y,z,e);
174  const double w = cst ? cst_val_2 : C(2,c,x,y,z,e);
175 
176  if (k == 0)
177  {
178  const double A11 = (J22 * J33) - (J23 * J32);
179  const double A12 = (J32 * J13) - (J12 * J33);
180  const double A13 = (J12 * J23) - (J22 * J13);
181  QQQ(z,y,x) = A11*u + A12*v + A13*w;
182 
183  }
184 
185  if (k == 1)
186  {
187  const double A21 = (J31 * J23) - (J21 * J33);
188  const double A22 = (J11 * J33) - (J13 * J31);
189  const double A23 = (J21 * J13) - (J11 * J23);
190  QQQ(z,y,x) = A21*u + A22*v + A23*w;
191  }
192 
193  if (k == 2)
194  {
195  const double A31 = (J21 * J32) - (J31 * J22);
196  const double A32 = (J31 * J12) - (J11 * J32);
197  const double A33 = (J11 * J22) - (J12 * J21);
198  QQQ(z,y,x) = A31*u + A32*v + A33*w;
199  }
200 
201  QQQ(z,y,x) *= W(x,y,z);
202  }
203  }
204  MFEM_SYNC_THREAD;
205  }
206  MFEM_FOREACH_THREAD(qz,x,q)
207  {
208  MFEM_FOREACH_THREAD(qy,y,q)
209  {
210  for (int dx = 0; dx < d; ++dx) { r_u[dx] = 0.0; }
211  for (int qx = 0; qx < q; ++qx)
212  {
213  const double r_v = QQQ(qz,qy,qx);
214  for (int dx = 0; dx < d; ++dx)
215  {
216  r_u[dx] += (k == 0 ? Gt(qx,dx) : Bt(qx,dx)) * r_v;
217  }
218  }
219  for (int dx = 0; dx < d; ++dx) { QQD(qz,qy,dx) = r_u[dx]; }
220  }
221  }
222  MFEM_SYNC_THREAD;
223  MFEM_FOREACH_THREAD(qz,y,q)
224  {
225  MFEM_FOREACH_THREAD(dx,x,d)
226  {
227  for (int dy = 0; dy < d; ++dy) { r_u[dy] = 0.0; }
228  for (int qy = 0; qy < q; ++qy)
229  {
230  const double r_v = QQD(qz,qy,dx);
231  for (int dy = 0; dy < d; ++dy)
232  {
233  r_u[dy] += (k == 1 ? Gt(qy,dy) : Bt(qy,dy)) * r_v;
234  }
235  }
236  for (int dy = 0; dy < d; ++dy) { QDD(qz,dy,dx) = r_u[dy]; }
237  }
238  }
239  MFEM_SYNC_THREAD;
240  MFEM_FOREACH_THREAD(dy,y,d)
241  {
242  MFEM_FOREACH_THREAD(dx,x,d)
243  {
244  for (int dz = 0; dz < d; ++dz) { r_u[dz] = 0.0; }
245  for (int qz = 0; qz < q; ++qz)
246  {
247  const double r_v = QDD(qz,dy,dx);
248  for (int dz = 0; dz < d; ++dz)
249  {
250  r_u[dz] += (k == 2 ? Gt(qz,dz) : Bt(qz,dz)) * r_v;
251  }
252  }
253  for (int dz = 0; dz < d; ++dz) { Y(dx,dy,dz,c,e) += r_u[dz]; }
254  }
255  }
256  MFEM_SYNC_THREAD;
257  } // dim
258  } // vdim
259  });
260 }
261 
262 static void DLFGradAssemble(const FiniteElementSpace &fes,
263  const IntegrationRule *ir,
264  const Array<int> &markers,
265  const Vector &coeff,
266  Vector &y)
267 {
268  Mesh *mesh = fes.GetMesh();
269  const int dim = mesh->Dimension();
270  const FiniteElement &el = *fes.GetFE(0);
271  const MemoryType mt = Device::GetDeviceMemoryType();
272  const DofToQuad &maps = el.GetDofToQuad(*ir, DofToQuad::TENSOR);
273  const int d = maps.ndof, q = maps.nqpt;
274  constexpr int flags = GeometricFactors::JACOBIANS;
275  const GeometricFactors *geom = mesh->GetGeometricFactors(*ir, flags, mt);
276  decltype(&DLFGradAssemble2D<>) ker =
277  dim == 2 ? DLFGradAssemble2D<> : DLFGradAssemble3D<>;
278 
279  if (dim==2)
280  {
281  if (d==1 && q==1) { ker=DLFGradAssemble2D<1,1>; }
282  if (d==2 && q==2) { ker=DLFGradAssemble2D<2,2>; }
283  if (d==3 && q==3) { ker=DLFGradAssemble2D<3,3>; }
284  if (d==4 && q==4) { ker=DLFGradAssemble2D<4,4>; }
285  if (d==5 && q==5) { ker=DLFGradAssemble2D<5,5>; }
286  if (d==2 && q==3) { ker=DLFGradAssemble2D<2,3>; }
287  if (d==3 && q==4) { ker=DLFGradAssemble2D<3,4>; }
288  if (d==4 && q==5) { ker=DLFGradAssemble2D<4,5>; }
289  if (d==5 && q==6) { ker=DLFGradAssemble2D<5,6>; }
290  }
291 
292  if (dim==3)
293  {
294  if (d==1 && q==1) { ker=DLFGradAssemble3D<1,1>; }
295  if (d==2 && q==2) { ker=DLFGradAssemble3D<2,2>; }
296  if (d==3 && q==3) { ker=DLFGradAssemble3D<3,3>; }
297  if (d==4 && q==4) { ker=DLFGradAssemble3D<4,4>; }
298  if (d==5 && q==5) { ker=DLFGradAssemble3D<5,5>; }
299  if (d==2 && q==3) { ker=DLFGradAssemble3D<2,3>; }
300  if (d==3 && q==4) { ker=DLFGradAssemble3D<3,4>; }
301  if (d==4 && q==5) { ker=DLFGradAssemble3D<4,5>; }
302  if (d==5 && q==6) { ker=DLFGradAssemble3D<5,6>; }
303  }
304 
305  MFEM_VERIFY(ker, "No kernel ndof " << d << " nqpt " << q);
306 
307  const int vdim = fes.GetVDim();
308  const int ne = fes.GetMesh()->GetNE();
309  const int *M = markers.Read();
310  const double *B = maps.B.Read();
311  const double *G = maps.G.Read();
312  const double *J = geom->J.Read();
313  const double *W = ir->GetWeights().Read();
314  double *Y = y.ReadWrite();
315  ker(vdim, ne, d, q, M, B, G, J, W, coeff, Y);
316 }
317 
318 void DomainLFGradIntegrator::AssembleDevice(const FiniteElementSpace &fes,
319  const Array<int> &markers,
320  Vector &b)
321 {
322 
323  const FiniteElement &fe = *fes.GetFE(0);
324  const int qorder = 2 * fe.GetOrder();
325  const Geometry::Type gtype = fe.GetGeomType();
326  const IntegrationRule *ir = IntRule ? IntRule : &IntRules.Get(gtype, qorder);
327 
328  QuadratureSpace qs(*fes.GetMesh(), *ir);
329  CoefficientVector coeff(Q, qs, CoefficientStorage::COMPRESSED);
330  DLFGradAssemble(fes, ir, markers, coeff, b);
331 }
332 
333 void VectorDomainLFGradIntegrator::AssembleDevice(const FiniteElementSpace &fes,
334  const Array<int> &markers,
335  Vector &b)
336 {
337  const FiniteElement &fe = *fes.GetFE(0);
338  const int qorder = 2 * fe.GetOrder();
339  const Geometry::Type gtype = fe.GetGeomType();
340  const IntegrationRule *ir = IntRule ? IntRule : &IntRules.Get(gtype, qorder);
341 
342  QuadratureSpace qs(*fes.GetMesh(), *ir);
343  CoefficientVector coeff(Q, qs, CoefficientStorage::COMPRESSED);
344  DLFGradAssemble(fes, ir, markers, coeff, b);
345 }
346 
347 } // namespace mfem
mfem::FiniteElement
Abstract class for all finite elements.
Definition: fe_base.hpp:235
mfem::IntegrationRule
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
mfem::IntegrationRules::Get
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:923
mfem::DofToQuad::TENSOR
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition: fe_base.hpp:162
mfem::FiniteElement::GetOrder
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe_base.hpp:327
mfem::LinearFormIntegrator::IntRule
const IntegrationRule * IntRule
Definition: lininteg.hpp:27
mfem::CoefficientVector
Class to represent a coefficient evaluated at quadrature points.
Definition: coefficient.hpp:2185
mfem::DeviceMatrix
DeviceTensor< 2, double > DeviceMatrix
Definition: dtensor.hpp:143
mfem::FiniteElement::GetGeomType
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:320
mfem::MAX_Q1D
const int MAX_Q1D
Definition: forall.hpp:29
mfem::FiniteElementSpace::GetMesh
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:441
b
double b
Definition: lissajous.cpp:42
mfem::Device::GetDeviceMemoryType
static MemoryType GetDeviceMemoryType()
Get the current Device MemoryType. This is the MemoryType used by most MFEM classes when allocating m...
Definition: device.hpp:273
mfem::CoefficientStorage::COMPRESSED
Enable all above compressions.
mfem::DeviceCube
DeviceTensor< 3, double > DeviceCube
Definition: dtensor.hpp:146
mfem::FiniteElementSpace
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
mfem::MemoryType
MemoryType
Memory types supported by MFEM.
Definition: mem_manager.hpp:31
mfem::VectorDomainLFGradIntegrator::AssembleDevice
virtual void AssembleDevice(const FiniteElementSpace &fes, const Array< int > &markers, Vector &b) override
Method defining assembly on device.
Definition: lininteg_domain_grad.cpp:333
mfem::DomainLFGradIntegrator::AssembleDevice
virtual void AssembleDevice(const FiniteElementSpace &fes, const Array< int > &markers, Vector &b)
Method defining assembly on device.
Definition: lininteg_domain_grad.cpp:318
dim
int dim
Definition: ex24.cpp:53
mfem::MAX_D1D
const int MAX_D1D
Definition: forall.hpp:28
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:2781
mfem::Vector
Vector data type.
Definition: vector.hpp:60
mfem::QuadratureSpace
Class representing the storage layout of a QuadratureFunction.
Definition: qspace.hpp:92
mfem::u
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
mfem::IntRules
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
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
Generated on Wed Dec 7 2022 10:41:46 for MFEM by   doxygen 1.8.5