MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
lininteg_domain_grad.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, 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 "../../fem/kernels.hpp"
14#include "../fem.hpp"
15
16namespace mfem
17{
18
19template<int T_D1D = 0, int T_Q1D = 0> static
20void DLFGradAssemble2D(const int vdim, const int ne, const int d, const int q,
21 const int *markers, const real_t *b, const real_t *g,
22 const real_t *jacobians,
23 const real_t *weights, const Vector &coeff, real_t *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(ne, q, q, [=] MFEM_HOST_DEVICE (int e)
36 {
37 if (M(e) == 0) { return; } // ignore
38
39 constexpr int Q = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
40 constexpr int D = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
41
42 MFEM_SHARED real_t sBGt[2][Q*D];
43 MFEM_SHARED real_t sQQ[2][Q*Q];
44 MFEM_SHARED real_t 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 real_t cst_val0 = C(0,c,0,0,0);
59 const real_t 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 real_t w = W(x,y);
66 const real_t J11 = J(x,y,0,0,e);
67 const real_t J21 = J(x,y,1,0,e);
68 const real_t J12 = J(x,y,0,1,e);
69 const real_t J22 = J(x,y,1,1,e);
70 const real_t u = cst ? cst_val0 : C(0,c,x,y,e);
71 const real_t 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 real_t 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 real_t 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
111template<int T_D1D = 0, int T_Q1D = 0> static
112void DLFGradAssemble3D(const int vdim, const int ne, const int d, const int q,
113 const int *markers, const real_t *b, const real_t *g,
114 const real_t *jacobians,
115 const real_t *weights, const Vector &coeff,
116 real_t *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(ne, q, q, [=] MFEM_HOST_DEVICE (int e)
130 {
131 if (M(e) == 0) { return; } // ignore
132
133 constexpr int Q = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
134 constexpr int D = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
135 constexpr int MQD = (Q >= D) ? Q : D;
136
137 MFEM_SHARED real_t sBGt[2][Q*D];
138 const DeviceMatrix Bt(sBGt[0], q,d), Gt(sBGt[1], q,d);
139
140 MFEM_SHARED real_t sQQQ[MQD*MQD*MQD];
141 const DeviceCube QQQ(sQQQ, MQD,MQD,MQD);
142
143 kernels::internal::LoadBGt<D,Q>(d,q,B,G,sBGt);
144
145 for (int c = 0; c < vdim; ++c)
146 {
147 const real_t cst_val_0 = C(0,c,0,0,0,0);
148 const real_t cst_val_1 = C(1,c,0,0,0,0);
149 const real_t cst_val_2 = C(2,c,0,0,0,0);
150
151 for (int k = 0; k < 3; ++k)
152 {
153 for (int z = 0; z < q; ++z)
154 {
155 MFEM_FOREACH_THREAD(y,y,q)
156 {
157 MFEM_FOREACH_THREAD(x,x,q)
158 {
159 const real_t J11 = J(x,y,z,0,0,e);
160 const real_t J21 = J(x,y,z,1,0,e);
161 const real_t J31 = J(x,y,z,2,0,e);
162 const real_t J12 = J(x,y,z,0,1,e);
163 const real_t J22 = J(x,y,z,1,1,e);
164 const real_t J32 = J(x,y,z,2,1,e);
165 const real_t J13 = J(x,y,z,0,2,e);
166 const real_t J23 = J(x,y,z,1,2,e);
167 const real_t J33 = J(x,y,z,2,2,e);
168
169 const real_t u = cst ? cst_val_0 : C(0,c,x,y,z,e);
170 const real_t v = cst ? cst_val_1 : C(1,c,x,y,z,e);
171 const real_t w = cst ? cst_val_2 : C(2,c,x,y,z,e);
172
173 if (k == 0)
174 {
175 const real_t A11 = (J22 * J33) - (J23 * J32);
176 const real_t A12 = (J32 * J13) - (J12 * J33);
177 const real_t A13 = (J12 * J23) - (J22 * J13);
178 QQQ(z,y,x) = A11*u + A12*v + A13*w;
179
180 }
181
182 if (k == 1)
183 {
184 const real_t A21 = (J31 * J23) - (J21 * J33);
185 const real_t A22 = (J11 * J33) - (J13 * J31);
186 const real_t A23 = (J21 * J13) - (J11 * J23);
187 QQQ(z,y,x) = A21*u + A22*v + A23*w;
188 }
189
190 if (k == 2)
191 {
192 const real_t A31 = (J21 * J32) - (J31 * J22);
193 const real_t A32 = (J31 * J12) - (J11 * J32);
194 const real_t A33 = (J11 * J22) - (J12 * J21);
195 QQQ(z,y,x) = A31*u + A32*v + A33*w;
196 }
197
198 QQQ(z,y,x) *= W(x,y,z);
199 }
200 }
201 MFEM_SYNC_THREAD;
202 }
203 MFEM_FOREACH_THREAD(qz,x,q)
204 {
205 MFEM_FOREACH_THREAD(qy,y,q)
206 {
207 real_t r_u[Q];
208 for (int qx = 0; qx < q; ++qx) { r_u[qx] = QQQ(qz,qy,qx); }
209 for (int dx = 0; dx < d; ++dx)
210 {
211 real_t u = 0.0;
212 for (int qx = 0; qx < q; ++qx)
213 {
214 u += (k == 0 ? Gt(qx,dx) : Bt(qx,dx)) * r_u[qx];
215 }
216 QQQ(qz,qy,dx) = u;
217 }
218 }
219 }
220 MFEM_SYNC_THREAD;
221 MFEM_FOREACH_THREAD(qz,y,q)
222 {
223 MFEM_FOREACH_THREAD(dx,x,d)
224 {
225 real_t r_u[Q];
226 for (int qy = 0; qy < q; ++qy) { r_u[qy] = QQQ(qz,qy,dx); }
227 for (int dy = 0; dy < d; ++dy)
228 {
229 real_t u = 0.0;
230 for (int qy = 0; qy < q; ++qy)
231 {
232 u += (k == 1 ? Gt(qy,dy) : Bt(qy,dy)) * r_u[qy];
233 }
234 QQQ(qz,dy,dx) = u;
235 }
236 }
237 }
238 MFEM_SYNC_THREAD;
239 MFEM_FOREACH_THREAD(dy,y,d)
240 {
241 MFEM_FOREACH_THREAD(dx,x,d)
242 {
243 real_t r_u[Q];
244 for (int qz = 0; qz < q; ++qz) { r_u[qz] = QQQ(qz,dy,dx); }
245 for (int dz = 0; dz < d; ++dz)
246 {
247 real_t u = 0.0;
248 for (int qz = 0; qz < q; ++qz)
249 {
250 u += (k == 2 ? Gt(qz,dz) : Bt(qz,dz)) * r_u[qz];
251 }
252 Y(dx,dy,dz,c,e) += u;
253 }
254 }
255 }
256 MFEM_SYNC_THREAD;
257 } // dim
258 } // vdim
259 });
260}
261
262static 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.GetTypicalFE();
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 real_t *B = maps.B.Read();
311 const real_t *G = maps.G.Read();
312 const real_t *J = geom->J.Read();
313 const real_t *W = ir->GetWeights().Read();
314 real_t *Y = y.ReadWrite();
315 ker(vdim, ne, d, q, M, B, G, J, W, coeff, Y);
316}
317
318void DomainLFGradIntegrator::AssembleDevice(const FiniteElementSpace &fes,
319 const Array<int> &markers,
320 Vector &b)
321{
322
323 const FiniteElement &fe = *fes.GetTypicalFE();
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
333void VectorDomainLFGradIntegrator::AssembleDevice(const FiniteElementSpace &fes,
334 const Array<int> &markers,
335 Vector &b)
336{
337 const FiniteElement &fe = *fes.GetTypicalFE();
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::CoefficientVector
Class to represent a coefficient evaluated at quadrature points.
Definition coefficient.hpp:2429
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:274
mfem::DofToQuad::TENSOR
@ TENSOR
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition fe_base.hpp:165
mfem::DomainLFGradIntegrator::AssembleDevice
void AssembleDevice(const FiniteElementSpace &fes, const Array< int > &markers, Vector &b) override
Method defining assembly on device.
Definition lininteg_domain_grad.cpp:318
mfem::FiniteElementSpace
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:244
mfem::FiniteElementSpace::GetMesh
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:679
mfem::FiniteElementSpace::GetTypicalFE
const FiniteElement * GetTypicalFE() const
Return GetFE(0) if the local mesh is not empty; otherwise return a typical FE based on the Geometry t...
Definition fespace.cpp:3871
mfem::FiniteElement
Abstract class for all finite elements.
Definition fe_base.hpp:244
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:338
mfem::FiniteElement::GetGeomType
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition fe_base.hpp:331
mfem::IntegrationRule
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
mfem::IntegrationRules::Get
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition intrules.cpp:1006
mfem::Integrator::IntRule
const IntegrationRule * IntRule
Definition integrator.hpp:70
mfem::QuadratureSpace
Class representing the storage layout of a QuadratureFunction.
Definition qspace.hpp:120
mfem::VectorDomainLFGradIntegrator::AssembleDevice
void AssembleDevice(const FiniteElementSpace &fes, const Array< int > &markers, Vector &b) override
Method defining assembly on device.
Definition lininteg_domain_grad.cpp:333
mfem::Vector
Vector data type.
Definition vector.hpp:82
dim
int dim
Definition ex24.cpp:53
b
real_t b
Definition lissajous.cpp:42
mfem::u
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
mfem::DeviceMatrix
DeviceTensor< 2, real_t > DeviceMatrix
Definition dtensor.hpp:143
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
mfem::CoefficientStorage::COMPRESSED
@ COMPRESSED
Enable all above compressions.
mfem::forall_2D
void forall_2D(int N, int X, int Y, lambda &&body)
Definition forall.hpp:762
mfem::real_t
float real_t
Definition config.hpp:43
mfem::MemoryType
MemoryType
Memory types supported by MFEM.
Definition mem_manager.hpp:39
mfem::DeviceCube
DeviceTensor< 3, real_t > DeviceCube
Definition dtensor.hpp:146
mfem::IntRules
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition intrules.hpp:492
Generated on Wed Apr 9 2025 16:53:04 for MFEM by doxygen 1.11.0