MFEM v4.9.0
Finite element discretization library
Loading...
Searching...
No Matches
bilininteg_vecmass_pa.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 "../bilininteg.hpp"
15
16#include "./bilininteg_vecmass_pa.hpp" // IWYU pragma: keep
17
18namespace mfem
19{
20
22{
23 Mesh *mesh = fes.GetMesh();
24 const FiniteElement &el = *fes.GetTypicalFE();
26 const auto *ir = IntRule ? IntRule : &MassIntegrator::GetRule(el, el, Trans);
27
28 if (DeviceCanUseCeed())
29 {
30 delete ceedOp;
31 const bool mixed =
32 mesh->GetNumGeometries(mesh->Dimension()) > 1 || fes.IsVariableOrder();
33 if (mixed) { ceedOp = new ceed::MixedPAMassIntegrator(*this, fes, Q); }
34 else { ceedOp = new ceed::PAMassIntegrator(fes, *ir, Q); }
35 return;
36 }
37
38 // If vdim is not set, set it to the space dimension
39 vdim = (vdim == -1) ? Trans.GetSpaceDim() : vdim;
40 MFEM_VERIFY(vdim == fes.GetVDim(), "vdim != fes.GetVDim()");
41 MFEM_VERIFY(vdim == mesh->Dimension(), "vdim != dim");
42
45 : pa_mt;
46
47 ne = mesh->GetNE();
48 dim = mesh->Dimension();
49 const int nq = ir->GetNPoints();
50 const int sdim = mesh->SpaceDimension();
53 dofs1D = maps->ndof;
54 quad1D = maps->nqpt;
55 const int q1d = quad1D;
56
57 if (!(dim == 2 || dim == 3)) { MFEM_ABORT("Dimension not supported."); }
58
59 QuadratureSpace qs(*mesh, *ir);
60 CoefficientVector coeff(qs);
61
62 if (Q)
63 {
64 coeff.Project(*Q);
65 }
66 else if (VQ)
67 {
68 coeff.Project(*VQ);
69 MFEM_VERIFY(VQ->GetVDim() == vdim, "VQ vdim vs. vdim error");
70 }
71 else if (MQ)
72 {
73 coeff.ProjectTranspose(*MQ);
74 MFEM_VERIFY(MQ->GetVDim() == vdim, "MQ dimension vs. vdim error");
75 MFEM_VERIFY(coeff.Size() == (vdim*vdim) * ne * nq, "MQ size error");
76 }
77 else { coeff.SetConstant(1.0); }
78
79 coeff_vdim = coeff.GetVDim();
80 const bool const_coeff = coeff_vdim == 1;
81 const bool vector_coeff = coeff_vdim == vdim;
82 const bool matrix_coeff = coeff_vdim == vdim * vdim;
83 MFEM_VERIFY(const_coeff + vector_coeff + matrix_coeff == 1, "");
84
85 pa_data.SetSize(coeff_vdim * nq * ne, mt);
86
87 const auto w_r = ir->GetWeights().Read();
88
89 if (dim == 2)
90 {
91 const auto W = Reshape(w_r, q1d, q1d);
92 const auto C = Reshape(coeff.Read(), coeff_vdim, q1d, q1d, ne);
93 const auto J = Reshape(geom->J.Read(), q1d, q1d, sdim, dim, ne);
94 auto D = Reshape(pa_data.Write(), q1d, q1d, coeff_vdim, ne);
95
96 mfem::forall_2D(ne, q1d, q1d, [=] MFEM_HOST_DEVICE(int e)
97 {
98 MFEM_FOREACH_THREAD(qy, y, q1d)
99 {
100 MFEM_FOREACH_THREAD(qx, x, q1d)
101 {
102 const real_t J11 = J(qx, qy, 0, 0, e), J12 = J(qx, qy, 1, 0, e);
103 const real_t J21 = J(qx, qy, 0, 1, e), J22 = J(qx, qy, 1, 1, e);
104 const real_t detJ = (J11 * J22) - (J21 * J12);
105 const real_t w_det = W(qx, qy) * detJ;
106 D(qx, qy, 0, e) = C(0, qx, qy, e) * w_det;
107 if (const_coeff) { continue; }
108 D(qx, qy, 1, e) = C(1, qx, qy, e) * w_det;
109 if (vector_coeff) { continue; }
110 assert(matrix_coeff);
111 D(qx, qy, 2, e) = C(2, qx, qy, e) * w_det;
112 D(qx, qy, 3, e) = C(3, qx, qy, e) * w_det;
113 }
114 }
115 });
116 }
117 else if (dim == 3)
118 {
119 const auto W = Reshape(w_r, q1d, q1d, q1d);
120 const auto C = Reshape(coeff.Read(), coeff_vdim, q1d, q1d, q1d, ne);
121 const auto J = Reshape(geom->J.Read(), q1d, q1d, q1d, sdim, dim, ne);
122 auto D = Reshape(pa_data.Write(), q1d, q1d, q1d, coeff_vdim, ne);
123
124 mfem::forall_3D(ne, q1d, q1d, q1d, [=] MFEM_HOST_DEVICE(int e)
125 {
126 MFEM_FOREACH_THREAD(qz, z, q1d)
127 {
128 MFEM_FOREACH_THREAD(qy, y, q1d)
129 {
130 MFEM_FOREACH_THREAD(qx, x, q1d)
131 {
132 const real_t J11 = J(qx, qy, qz, 0, 0, e),
133 J12 = J(qx, qy, qz, 0, 1, e),
134 J13 = J(qx, qy, qz, 0, 2, e);
135 const real_t J21 = J(qx, qy, qz, 1, 0, e),
136 J22 = J(qx, qy, qz, 1, 1, e),
137 J23 = J(qx, qy, qz, 1, 2, e);
138 const real_t J31 = J(qx, qy, qz, 2, 0, e),
139 J32 = J(qx, qy, qz, 2, 1, e),
140 J33 = J(qx, qy, qz, 2, 2, e);
141 const real_t detJ = J11 * (J22 * J33 - J32 * J23) -
142 J21 * (J12 * J33 - J32 * J13) +
143 J31 * (J12 * J23 - J22 * J13);
144 const real_t w_det = W(qx, qy, qz) * detJ;
145 D(qx, qy, qz, 0, e) = C(0, qx, qy, qz, e) * w_det;
146 if (const_coeff) { continue; }
147 D(qx, qy, qz, 1, e) = C(1, qx, qy, qz, e) * w_det;
148 D(qx, qy, qz, 2, e) = C(2, qx, qy, qz, e) * w_det;
149 if (vector_coeff) { continue; }
150 D(qx, qy, qz, 3, e) = C(3, qx, qy, qz, e) * w_det;
151 D(qx, qy, qz, 4, e) = C(4, qx, qy, qz, e) * w_det;
152 D(qx, qy, qz, 5, e) = C(5, qx, qy, qz, e) * w_det;
153 D(qx, qy, qz, 6, e) = C(6, qx, qy, qz, e) * w_det;
154 D(qx, qy, qz, 7, e) = C(7, qx, qy, qz, e) * w_det;
155 D(qx, qy, qz, 8, e) = C(8, qx, qy, qz, e) * w_det;
156 }
157 }
158 }
159 });
160 }
161 else
162 {
163 MFEM_ABORT("Unknown VectorMassIntegrator::AssemblePA kernel for"
164 << " dim:" << dim << ", vdim:" << vdim << ", sdim:" << sdim);
165 }
166}
167
169{
170 // Use CEED backend if available
171 if (DeviceCanUseCeed()) { return ceedOp->AddMult(x, y); }
172
173 // Add the VectorMassAddMultPA specializations
174 static const auto vector_mass_kernel_specializations =
175 ( // 2D
176 VectorMassIntegrator::VectorMassAddMultPA::Specialization<2, 2,2>::Add(),
177 VectorMassIntegrator::VectorMassAddMultPA::Specialization<2, 3,3>::Add(),
178 VectorMassIntegrator::VectorMassAddMultPA::Specialization<2, 3,4>::Add(),
179 VectorMassIntegrator::VectorMassAddMultPA::Specialization<2, 4,4>::Add(),
180 VectorMassIntegrator::VectorMassAddMultPA::Specialization<2, 4,6>::Add(),
181 VectorMassIntegrator::VectorMassAddMultPA::Specialization<2, 5,5>::Add(),
182 VectorMassIntegrator::VectorMassAddMultPA::Specialization<2, 6,6>::Add(),
183 VectorMassIntegrator::VectorMassAddMultPA::Specialization<2, 7,7>::Add(),
184 VectorMassIntegrator::VectorMassAddMultPA::Specialization<2, 8,8>::Add(),
185 VectorMassIntegrator::VectorMassAddMultPA::Specialization<2, 9,9>::Add(),
186 // 3D
187 VectorMassIntegrator::VectorMassAddMultPA::Specialization<3, 2,2>::Add(),
188 VectorMassIntegrator::VectorMassAddMultPA::Specialization<3, 2,3>::Add(),
189 VectorMassIntegrator::VectorMassAddMultPA::Specialization<3, 3,4>::Add(),
190 VectorMassIntegrator::VectorMassAddMultPA::Specialization<3, 3,5>::Add(),
191 VectorMassIntegrator::VectorMassAddMultPA::Specialization<3, 4,5>::Add(),
192 VectorMassIntegrator::VectorMassAddMultPA::Specialization<3, 4,6>::Add(),
193 VectorMassIntegrator::VectorMassAddMultPA::Specialization<3, 4,8>::Add(),
194 VectorMassIntegrator::VectorMassAddMultPA::Specialization<3, 5,6>::Add(),
195 VectorMassIntegrator::VectorMassAddMultPA::Specialization<3, 5,8>::Add(),
196 VectorMassIntegrator::VectorMassAddMultPA::Specialization<3, 6,7>::Add(),
197 VectorMassIntegrator::VectorMassAddMultPA::Specialization<3, 7,8>::Add(),
198 VectorMassIntegrator::VectorMassAddMultPA::Specialization<3, 8,9>::Add(),
199 true);
200 MFEM_CONTRACT_VAR(vector_mass_kernel_specializations);
201
202 VectorMassAddMultPA::Run(dim, dofs1D, quad1D,
203 ne, coeff_vdim, maps->B, pa_data, x, y,
204 dofs1D, quad1D);
205
206}
207
208template <const int T_D1D = 0, const int T_Q1D = 0>
209static void PAVectorMassAssembleDiagonal2D(const int NE,
210 const Array<real_t> &b,
211 const Vector &pa_data, Vector &diag,
212 const int d1d = 0, const int q1d = 0)
213{
214 constexpr int VDIM = 2;
215 const int D1D = T_D1D ? T_D1D : d1d;
216 const int Q1D = T_Q1D ? T_Q1D : q1d;
217 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
218 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
219 const auto B = Reshape(b.Read(), Q1D, D1D);
220 const auto D = Reshape(pa_data.Read(), Q1D, Q1D, NE);
221 auto Y = Reshape(diag.ReadWrite(), D1D, D1D, VDIM, NE);
222
223 mfem::forall(NE, [=] MFEM_HOST_DEVICE(int e)
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 : DofQuadLimits::MAX_D1D;
228 constexpr int max_Q1D = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
229
230 real_t temp[max_Q1D][max_D1D];
231 for (int qx = 0; qx < Q1D; ++qx)
232 {
233 for (int dy = 0; dy < D1D; ++dy)
234 {
235 temp[qx][dy] = 0.0;
236 for (int qy = 0; qy < Q1D; ++qy)
237 {
238 temp[qx][dy] += B(qy, dy) * B(qy, dy) * D(qx, qy, e);
239 }
240 }
241 }
242 for (int dy = 0; dy < D1D; ++dy)
243 {
244 for (int dx = 0; dx < D1D; ++dx)
245 {
246 real_t temp1 = 0.0;
247 for (int qx = 0; qx < Q1D; ++qx)
248 {
249 temp1 += B(qx, dx) * B(qx, dx) * temp[qx][dy];
250 }
251 Y(dx, dy, 0, e) = temp1;
252 Y(dx, dy, 1, e) = temp1;
253 }
254 }
255 });
256}
257
258template <const int T_D1D = 0, const int T_Q1D = 0>
259static void PAVectorMassAssembleDiagonal3D(const int NE,
260 const Array<real_t> &B_,
261 const Vector &pa_data, Vector &diag,
262 const int d1d = 0, const int q1d = 0)
263{
264 constexpr int VDIM = 3;
265 const int D1D = T_D1D ? T_D1D : d1d;
266 const int Q1D = T_Q1D ? T_Q1D : q1d;
267 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
268 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
269 const auto B = Reshape(B_.Read(), Q1D, D1D);
270 MFEM_VERIFY(pa_data.Size() == Q1D * Q1D * Q1D * NE, "pa_data size error");
271 const auto D = Reshape(pa_data.Read(), Q1D, Q1D, Q1D, NE);
272 auto Y = Reshape(diag.ReadWrite(), D1D, D1D, D1D, VDIM, NE);
273 mfem::forall(NE, [=] MFEM_HOST_DEVICE(int e)
274 {
275 const int D1D = T_D1D ? T_D1D : d1d;
276 const int Q1D = T_Q1D ? T_Q1D : q1d;
277 // the following variables are evaluated at compile time
278 constexpr int max_D1D = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
279 constexpr int max_Q1D = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
280
281 real_t temp[max_Q1D][max_Q1D][max_D1D];
282 for (int qx = 0; qx < Q1D; ++qx)
283 {
284 for (int qy = 0; qy < Q1D; ++qy)
285 {
286 for (int dz = 0; dz < D1D; ++dz)
287 {
288 temp[qx][qy][dz] = 0.0;
289 for (int qz = 0; qz < Q1D; ++qz)
290 {
291 temp[qx][qy][dz] +=
292 B(qz, dz) * B(qz, dz) * D(qx, qy, qz, e);
293 }
294 }
295 }
296 }
297 real_t temp2[max_Q1D][max_D1D][max_D1D];
298 for (int qx = 0; qx < Q1D; ++qx)
299 {
300 for (int dz = 0; dz < D1D; ++dz)
301 {
302 for (int dy = 0; dy < D1D; ++dy)
303 {
304 temp2[qx][dy][dz] = 0.0;
305 for (int qy = 0; qy < Q1D; ++qy)
306 {
307 temp2[qx][dy][dz] +=
308 B(qy, dy) * B(qy, dy) * temp[qx][qy][dz];
309 }
310 }
311 }
312 }
313 for (int dz = 0; dz < D1D; ++dz)
314 {
315 for (int dy = 0; dy < D1D; ++dy)
316 {
317 for (int dx = 0; dx < D1D; ++dx)
318 {
319 real_t temp3 = 0.0;
320 for (int qx = 0; qx < Q1D; ++qx)
321 {
322 temp3 += B(qx, dx) * B(qx, dx) * temp2[qx][dy][dz];
323 }
324 Y(dx, dy, dz, 0, e) = temp3;
325 Y(dx, dy, dz, 1, e) = temp3;
326 Y(dx, dy, dz, 2, e) = temp3;
327 }
328 }
329 }
330 });
331}
332
333static void PAVectorMassAssembleDiagonal(const int dim, const int D1D,
334 const int Q1D, const int NE,
335 const Array<real_t> &B,
336 const Vector &pa_data,
337 Vector &diag)
338{
339 if (dim == 2)
340 {
341 return PAVectorMassAssembleDiagonal2D(NE, B, pa_data, diag, D1D, Q1D);
342 }
343 else if (dim == 3)
344 {
345 return PAVectorMassAssembleDiagonal3D(NE, B, pa_data, diag, D1D, Q1D);
346 }
347 MFEM_ABORT("Dimension not implemented.");
348}
349
351{
352 if (DeviceCanUseCeed()) { ceedOp->GetDiagonal(diag); }
353 else
354 {
355 MFEM_VERIFY(coeff_vdim == 1, "coeff_vdim != 1");
356 MFEM_VERIFY(!VQ && !MQ, "VQ and MQ not supported");
357 PAVectorMassAssembleDiagonal(dim, dofs1D, quad1D, ne, maps->B, pa_data, diag);
358 }
359}
360
361} // namespace mfem
Class to represent a coefficient evaluated at quadrature points.
void SetConstant(real_t constant)
Set this vector to the given constant.
void Project(Coefficient &coeff)
Evaluate the given Coefficient at the quadrature points defined by qs.
int GetVDim() const
Return the number of values per quadrature point.
void ProjectTranspose(MatrixCoefficient &coeff)
Project the transpose of coeff.
static MemoryType GetDeviceMemoryType()
Get the current Device MemoryType. This is the MemoryType used by most MFEM classes when allocating m...
Definition device.hpp:277
@ TENSOR
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition fe_base.hpp:165
Array< real_t > B
Basis functions evaluated at quadrature points.
Definition fe_base.hpp:193
int ndof
Number of degrees of freedom = number of basis functions. When mode is TENSOR, this is the 1D number.
Definition fe_base.hpp:178
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Definition fe_base.hpp:182
virtual int GetSpaceDim() const =0
Get the dimension of the target (physical) space.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:208
bool IsVariableOrder() const
Returns true if the space contains elements of varying polynomial orders.
Definition fespace.hpp:678
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:644
int GetVDim() const
Returns the vector dimension of the finite element space.
Definition fespace.hpp:819
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:3860
Abstract class for all finite elements.
Definition fe_base.hpp:247
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_base.cpp:375
Vector J
Jacobians of the element transformations at all quadrature points.
Definition mesh.hpp:3115
const IntegrationRule * IntRule
static const IntegrationRule & GetRule(const FiniteElement &trial_fe, const FiniteElement &test_fe, const ElementTransformation &Trans)
int GetVDim() const
For backward compatibility get the width of the matrix.
Mesh data type.
Definition mesh.hpp:65
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1377
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1306
ElementTransformation * GetTypicalElementTransformation()
If the local mesh is not empty return GetElementTransformation(0); otherwise, return the identity tra...
Definition mesh.cpp:393
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition mesh.hpp:1309
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:883
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
Definition mesh.cpp:7558
Class representing the storage layout of a QuadratureFunction.
Definition qspace.hpp:164
int GetVDim()
Returns dimension of the vector.
VectorCoefficient * VQ
void AssembleDiagonalPA(Vector &diag) override
Assemble diagonal and add it to Vector diag.
const DofToQuad * maps
Not owned.
void AddMultPA(const Vector &x, Vector &y) const override
Method for partially assembled action.
MatrixCoefficient * MQ
void AssemblePA(const FiniteElementSpace &fes) override
Method defining partial assembly.
const GeometricFactors * geom
Not owned.
Vector data type.
Definition vector.hpp:82
virtual const real_t * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:520
virtual real_t * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:536
int Size() const
Returns the size of the vector.
Definition vector.hpp:234
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:584
virtual real_t * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:528
void GetDiagonal(mfem::Vector &diag) const
Definition operator.cpp:104
void AddMult(const mfem::Vector &x, mfem::Vector &y, const real_t a=1.0) const override
Operator application: y+=A(x) (default) or y+=a*A(x).
Definition operator.cpp:72
Represent a MassIntegrator with AssemblyLevel::Partial using libCEED.
Definition mass.hpp:27
int dim
Definition ex24.cpp:53
real_t b
Definition lissajous.cpp:42
mfem::real_t real_t
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:138
bool DeviceCanUseCeed()
Function that determines if a CEED kernel should be used, based on the current mfem::Device configura...
Definition util.cpp:33
void forall_2D(int N, int X, int Y, lambda &&body)
Definition forall.hpp:925
void forall_3D(int N, int X, int Y, int Z, lambda &&body)
Definition forall.hpp:937
float real_t
Definition config.hpp:46
MemoryType
Memory types supported by MFEM.
void forall(int N, lambda &&body)
Definition forall.hpp:839
static const DeviceDofQuadLimits & Get()
Return a const reference to the DeviceDofQuadLimits singleton.
Definition forall.hpp:128