MFEM v4.8.0
Finite element discretization library
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
lininteg_domain_vectorfe.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>
20static void HdivDLFAssemble2D(
21 const int ne, const int d, const int q, const int *markers, const real_t *bo,
22 const real_t *bc, const real_t *j, const real_t *weights,
23 const Vector &coeff, real_t *y)
24{
25 MFEM_VERIFY(T_D1D || d <= DeviceDofQuadLimits::Get().HDIV_MAX_D1D,
26 "Problem size too large.");
27 MFEM_VERIFY(T_Q1D || q <= DeviceDofQuadLimits::Get().HDIV_MAX_Q1D,
28 "Problem size too large.");
29
30 static constexpr int vdim = 2;
31 const auto F = coeff.Read();
32 const auto M = Reshape(markers, ne);
33 const auto BO = Reshape(bo, q, d-1);
34 const auto BC = Reshape(bc, q, d);
35 const auto J = Reshape(j, q, q, vdim, vdim, ne);
36 const auto W = Reshape(weights, q, q);
37 const bool cst = coeff.Size() == vdim;
38 const auto C = cst ? Reshape(F,vdim,1,1,1) : Reshape(F,vdim,q,q,ne);
39 auto Y = Reshape(y, 2*(d-1)*d, ne);
40
41 mfem::forall_3D(ne, q, q, vdim, [=] MFEM_HOST_DEVICE (int e)
42 {
43 if (M(e) == 0) { return; } // ignore
44
45 constexpr int Q = T_Q1D ? T_Q1D : DofQuadLimits::HDIV_MAX_Q1D;
46 constexpr int D = T_D1D ? T_D1D : DofQuadLimits::HDIV_MAX_D1D;
47
48 MFEM_SHARED real_t sBot[Q*D];
49 MFEM_SHARED real_t sBct[Q*D];
50 MFEM_SHARED real_t sQQ[vdim*Q*Q];
51 MFEM_SHARED real_t sQD[vdim*Q*D];
52
53 // Bo and Bc into shared memory
54 const DeviceMatrix Bot(sBot, d-1, q);
55 kernels::internal::LoadB<D,Q>(d-1, q, BO, sBot);
56 const DeviceMatrix Bct(sBct, d, q);
57 kernels::internal::LoadB<D,Q>(d, q, BC, sBct);
58
59 const DeviceCube QQ(sQQ, q, q, vdim);
60 const DeviceCube QD(sQD, q, d, vdim);
61
62 MFEM_FOREACH_THREAD(vd,z,vdim)
63 {
64 const real_t cst_val_0 = C(0,0,0,0);
65 const real_t cst_val_1 = C(1,0,0,0);
66 MFEM_FOREACH_THREAD(y,y,q)
67 {
68 MFEM_FOREACH_THREAD(x,x,q)
69 {
70 const real_t J0 = J(x,y,0,vd,e);
71 const real_t J1 = J(x,y,1,vd,e);
72 const real_t C0 = cst ? cst_val_0 : C(0,x,y,e);
73 const real_t C1 = cst ? cst_val_1 : C(1,x,y,e);
74 QQ(x,y,vd) = W(x,y)*(J0*C0 + J1*C1);
75 }
76 }
77 }
78 MFEM_SYNC_THREAD;
79 MFEM_FOREACH_THREAD(vd,z,vdim)
80 {
81 const int nx = (vd == 0) ? d : d-1;
82 DeviceMatrix Btx = (vd == 0) ? Bct : Bot;
83 MFEM_FOREACH_THREAD(qy,y,q)
84 {
85 MFEM_FOREACH_THREAD(dx,x,nx)
86 {
87 real_t qd = 0.0;
88 for (int qx = 0; qx < q; ++qx)
89 {
90 qd += QQ(qx,qy,vd) * Btx(dx,qx);
91 }
92 QD(dx,qy,vd) = qd;
93 }
94 }
95 }
96 MFEM_SYNC_THREAD;
97 MFEM_FOREACH_THREAD(vd,z,vdim)
98 {
99 const int nx = (vd == 0) ? d : d-1;
100 const int ny = (vd == 1) ? d : d-1;
101 DeviceMatrix Bty = (vd == 1) ? Bct : Bot;
102 DeviceTensor<4> Yxy(Y, nx, ny, vdim, ne);
103 MFEM_FOREACH_THREAD(dy,y,ny)
104 {
105 MFEM_FOREACH_THREAD(dx,x,nx)
106 {
107 real_t dd = 0.0;
108 for (int qy = 0; qy < q; ++qy)
109 {
110 dd += QD(dx,qy,vd) * Bty(dy,qy);
111 }
112 Yxy(dx,dy,vd,e) += dd;
113 }
114 }
115 }
116 MFEM_SYNC_THREAD;
117 });
118}
119
120template<int T_D1D = 0, int T_Q1D = 0>
121static void HdivDLFAssemble3D(
122 const int ne, const int d, const int q, const int *markers, const real_t *bo,
123 const real_t *bc, const real_t *j, const real_t *weights,
124 const Vector &coeff, real_t *y)
125{
126 MFEM_VERIFY(T_D1D || d <= DeviceDofQuadLimits::Get().HDIV_MAX_D1D,
127 "Problem size too large.");
128 MFEM_VERIFY(T_Q1D || q <= DeviceDofQuadLimits::Get().HDIV_MAX_Q1D,
129 "Problem size too large.");
130
131 static constexpr int vdim = 3;
132 const auto F = coeff.Read();
133 const auto M = Reshape(markers, ne);
134 const auto BO = Reshape(bo, q, d-1);
135 const auto BC = Reshape(bc, q, d);
136 const auto J = Reshape(j, q, q, q, vdim, vdim, ne);
137 const auto W = Reshape(weights, q, q, q);
138 const bool cst = coeff.Size() == vdim;
139 const auto C = cst ? Reshape(F,vdim,1,1,1,1) : Reshape(F,vdim,q,q,q,ne);
140 auto Y = Reshape(y, 2*(d-1)*(d-1)*d, ne);
141
142 mfem::forall_3D(ne, q, q, vdim, [=] MFEM_HOST_DEVICE (int e)
143 {
144 if (M(e) == 0) { return; } // ignore
145
146 constexpr int Q = T_Q1D ? T_Q1D : DofQuadLimits::HDIV_MAX_Q1D;
147 constexpr int D = T_D1D ? T_D1D : DofQuadLimits::HDIV_MAX_D1D;
148
149 MFEM_SHARED real_t sBot[Q*D];
150 MFEM_SHARED real_t sBct[Q*D];
151
152 // Bo and Bc into shared memory
153 const DeviceMatrix Bot(sBot, d-1, q);
154 kernels::internal::LoadB<D,Q>(d-1, q, BO, sBot);
155 const DeviceMatrix Bct(sBct, d, q);
156 kernels::internal::LoadB<D,Q>(d, q, BC, sBct);
157
158 MFEM_SHARED real_t sm0[vdim*Q*Q*Q];
159 MFEM_SHARED real_t sm1[vdim*Q*Q*Q];
160 DeviceTensor<4> QQQ(sm1, q, q, q, vdim);
161 DeviceTensor<4> DQQ(sm0, d, q, q, vdim);
162 DeviceTensor<4> DDQ(sm1, d, d, q, vdim);
163
164 MFEM_FOREACH_THREAD(vd,z,vdim)
165 {
166 const real_t cst_val_0 = C(0,0,0,0,0);
167 const real_t cst_val_1 = C(1,0,0,0,0);
168 const real_t cst_val_2 = C(2,0,0,0,0);
169 MFEM_FOREACH_THREAD(y,y,q)
170 {
171 MFEM_FOREACH_THREAD(x,x,q)
172 {
173 for (int z = 0; z < q; ++z)
174 {
175 const real_t J0 = J(x,y,z,0,vd,e);
176 const real_t J1 = J(x,y,z,1,vd,e);
177 const real_t J2 = J(x,y,z,2,vd,e);
178 const real_t C0 = cst ? cst_val_0 : C(0,x,y,z,e);
179 const real_t C1 = cst ? cst_val_1 : C(1,x,y,z,e);
180 const real_t C2 = cst ? cst_val_2 : C(2,x,y,z,e);
181 QQQ(x,y,z,vd) = W(x,y,z)*(J0*C0 + J1*C1 + J2*C2);
182 }
183 }
184 }
185 }
186 MFEM_SYNC_THREAD;
187 // Apply Bt operator
188 MFEM_FOREACH_THREAD(vd,z,vdim)
189 {
190 const int nx = (vd == 0) ? d : d-1;
191 DeviceMatrix Btx = (vd == 0) ? Bct : Bot;
192 MFEM_FOREACH_THREAD(qy,y,q)
193 {
194 MFEM_FOREACH_THREAD(dx,x,nx)
195 {
196 real_t u[Q];
197 MFEM_UNROLL(Q)
198 for (int qz = 0; qz < q; ++qz) { u[qz] = 0.0; }
199 MFEM_UNROLL(Q)
200 for (int qx = 0; qx < q; ++qx)
201 {
202 MFEM_UNROLL(Q)
203 for (int qz = 0; qz < q; ++qz)
204 {
205 u[qz] += QQQ(qx,qy,qz,vd) * Btx(dx,qx);
206 }
207 }
208 MFEM_UNROLL(Q)
209 for (int qz = 0; qz < q; ++qz) { DQQ(dx,qy,qz,vd) = u[qz]; }
210 }
211 }
212 }
213 MFEM_SYNC_THREAD;
214 MFEM_FOREACH_THREAD(vd,z,vdim)
215 {
216 const int nx = (vd == 0) ? d : d-1;
217 const int ny = (vd == 1) ? d : d-1;
218 DeviceMatrix Bty = (vd == 1) ? Bct : Bot;
219 MFEM_FOREACH_THREAD(dy,y,ny)
220 {
221 MFEM_FOREACH_THREAD(dx,x,nx)
222 {
223 real_t u[Q];
224 MFEM_UNROLL(Q)
225 for (int qz = 0; qz < q; ++qz) { u[qz] = 0.0; }
226 MFEM_UNROLL(Q)
227 for (int qy = 0; qy < q; ++qy)
228 {
229 MFEM_UNROLL(Q)
230 for (int qz = 0; qz < q; ++qz)
231 {
232 u[qz] += DQQ(dx,qy,qz,vd) * Bty(dy,qy);
233 }
234 }
235 MFEM_UNROLL(Q)
236 for (int qz = 0; qz < q; ++qz) { DDQ(dx,dy,qz,vd) = u[qz]; }
237 }
238 }
239 }
240 MFEM_SYNC_THREAD;
241 MFEM_FOREACH_THREAD(vd,z,vdim)
242 {
243 const int nx = (vd == 0) ? d : d-1;
244 const int ny = (vd == 1) ? d : d-1;
245 const int nz = (vd == 2) ? d : d-1;
246 DeviceTensor<5> Yxyz(Y, nx, ny, nz, vdim, ne);
247 DeviceMatrix Btz = (vd == 2) ? Bct : Bot;
248 MFEM_FOREACH_THREAD(dy,y,ny)
249 {
250 MFEM_FOREACH_THREAD(dx,x,nx)
251 {
252 real_t u[D];
253 MFEM_UNROLL(D)
254 for (int dz = 0; dz < nz; ++dz) { u[dz] = 0.0; }
255 MFEM_UNROLL(Q)
256 for (int qz = 0; qz < q; ++qz)
257 {
258 MFEM_UNROLL(D)
259 for (int dz = 0; dz < nz; ++dz)
260 {
261 u[dz] += DDQ(dx,dy,qz,vd) * Btz(dz,qz);
262 }
263 }
264 MFEM_UNROLL(D)
265 for (int dz = 0; dz < nz; ++dz) { Yxyz(dx,dy,dz,vd,e) += u[dz]; }
266 }
267 }
268 }
269 MFEM_SYNC_THREAD;
270 });
271}
272
273static void HdivDLFAssemble(const FiniteElementSpace &fes,
274 const IntegrationRule *ir,
275 const Array<int> &markers,
276 const Vector &coeff,
277 Vector &y)
278{
279 Mesh &mesh = *fes.GetMesh();
280 const int dim = mesh.Dimension();
281 const FiniteElement *el = fes.GetTypicalFE();
282 const auto *vel = dynamic_cast<const VectorTensorFiniteElement *>(el);
283 MFEM_VERIFY(vel != nullptr, "Must be VectorTensorFiniteElement");
284 const MemoryType mt = Device::GetDeviceMemoryType();
285 const DofToQuad &maps_o = vel->GetDofToQuadOpen(*ir, DofToQuad::TENSOR);
286 const DofToQuad &maps_c = vel->GetDofToQuad(*ir, DofToQuad::TENSOR);
287 const int d = maps_c.ndof, q = maps_c.nqpt;
288 constexpr int flags = GeometricFactors::JACOBIANS;
289 const GeometricFactors *geom = mesh.GetGeometricFactors(*ir, flags, mt);
290 decltype(&HdivDLFAssemble2D<>) ker =
291 dim == 2 ? HdivDLFAssemble2D<> : HdivDLFAssemble3D<>;
292
293 if (dim==2)
294 {
295 if (d==1 && q==1) { ker=HdivDLFAssemble2D<1,1>; }
296 if (d==2 && q==2) { ker=HdivDLFAssemble2D<2,2>; }
297 if (d==3 && q==3) { ker=HdivDLFAssemble2D<3,3>; }
298 if (d==4 && q==4) { ker=HdivDLFAssemble2D<4,4>; }
299 if (d==5 && q==5) { ker=HdivDLFAssemble2D<5,5>; }
300 if (d==6 && q==6) { ker=HdivDLFAssemble2D<6,6>; }
301 if (d==7 && q==7) { ker=HdivDLFAssemble2D<7,7>; }
302 if (d==8 && q==8) { ker=HdivDLFAssemble2D<8,8>; }
303 }
304
305 if (dim==3)
306 {
307 if (d==2 && q==2) { ker=HdivDLFAssemble3D<2,2>; }
308 if (d==3 && q==3) { ker=HdivDLFAssemble3D<3,3>; }
309 if (d==4 && q==4) { ker=HdivDLFAssemble3D<4,4>; }
310 if (d==5 && q==5) { ker=HdivDLFAssemble3D<5,5>; }
311 if (d==6 && q==6) { ker=HdivDLFAssemble3D<6,6>; }
312 if (d==7 && q==7) { ker=HdivDLFAssemble3D<7,7>; }
313 if (d==8 && q==8) { ker=HdivDLFAssemble3D<8,8>; }
314 }
315
316 MFEM_VERIFY(ker, "No kernel ndof " << d << " nqpt " << q);
317
318 const int ne = mesh.GetNE();
319 const int *M = markers.Read();
320 const real_t *Bo = maps_o.B.Read();
321 const real_t *Bc = maps_c.B.Read();
322 const real_t *J = geom->J.Read();
323 const real_t *W = ir->GetWeights().Read();
324 real_t *Y = y.ReadWrite();
325 ker(ne, d, q, M, Bo, Bc, J, W, coeff, Y);
326}
327
328void VectorFEDomainLFIntegrator::AssembleDevice(const FiniteElementSpace &fes,
329 const Array<int> &markers,
330 Vector &b)
331{
332 const FiniteElement &fe = *fes.GetTypicalFE();
333 const int qorder = 2 * fe.GetOrder();
334 const Geometry::Type gtype = fe.GetGeomType();
335 const IntegrationRule *ir = IntRule ? IntRule : &IntRules.Get(gtype, qorder);
336
337 QuadratureSpace qs(*fes.GetMesh(), *ir);
338 CoefficientVector coeff(QF, qs, CoefficientStorage::COMPRESSED);
339
340 const int fe_type = fe.GetDerivType();
341 if (fe_type == FiniteElement::DIV)
342 {
343 HdivDLFAssemble(fes, ir, markers, coeff, b);
344 }
345 else
346 {
347 MFEM_ABORT("Not implemented.");
348 }
349}
350
351} // 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::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::GetDerivType
int GetDerivType() const
Returns the FiniteElement::DerivType of the element describing the spatial derivative method implemen...
Definition fe_base.hpp:365
mfem::FiniteElement::GetGeomType
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition fe_base.hpp:331
mfem::FiniteElement::DIV
@ DIV
Implements CalcDivShape methods.
Definition fe_base.hpp:306
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::VectorFEDomainLFIntegrator::AssembleDevice
void AssembleDevice(const FiniteElementSpace &fes, const Array< int > &markers, Vector &b) override
Method defining assembly on device.
Definition lininteg_domain_vectorfe.cpp:328
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_3D
void forall_3D(int N, int X, int Y, int Z, lambda &&body)
Definition forall.hpp:774
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
vel
void vel(const Vector &x, real_t t, Vector &u)
Definition navier_3dfoc.cpp:28
mfem::DeviceDofQuadLimits::Get
static const DeviceDofQuadLimits & Get()
Return a const reference to the DeviceDofQuadLimits singleton.
Definition forall.hpp:128
Generated on Wed Apr 9 2025 16:53:04 for MFEM by doxygen 1.11.0