MFEM v4.10.0
Finite element discretization library
Loading...
Searching...
No Matches
bilininteg_vecdiv_pa.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2026, 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 "./bilininteg_vecdiv_pa.hpp" // IWYU pragma: keep
15
16namespace mfem
17{
18
19// PA Divergence Assemble 2D kernel
20static void PADivergenceSetup2D(const int Q1D,
21 const int NE,
22 const Array<real_t> &w,
23 const Vector &j,
24 const real_t COEFF,
25 Vector &op)
26{
27 const auto W = Reshape(w.Read(), Q1D, Q1D);
28 const auto J = Reshape(j.Read(), Q1D, Q1D, 2, 2, NE);
29 auto y = Reshape(op.Write(), Q1D, Q1D, 2, 2, NE);
30
31 mfem::forall_2D(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
32 {
33 MFEM_FOREACH_THREAD_DIRECT(qx, x, Q1D)
34 {
35 MFEM_FOREACH_THREAD_DIRECT(qy, y, Q1D)
36 {
37 const real_t J11 = J(qx, qy, 0, 0, e);
38 const real_t J21 = J(qx, qy, 1, 0, e);
39 const real_t J12 = J(qx, qy, 0, 1, e);
40 const real_t J22 = J(qx, qy, 1, 1, e);
41 const real_t cw = W(qx, qy) * COEFF;
42 y(qx, qy, 0, 0, e) = cw * J22;
43 y(qx, qy, 0, 1, e) = cw * -J12;
44 y(qx, qy, 1, 0, e) = cw * -J21;
45 y(qx, qy, 1, 1, e) = cw * J11;
46 }
47 }
48 });
49}
50
51// PA Divergence Assemble 3D kernel
52static void PADivergenceSetup3D(const int Q1D,
53 const int NE,
54 const Array<real_t> &w,
55 const Vector &j,
56 const real_t COEFF,
57 Vector &op)
58{
59 const auto W = Reshape(w.Read(), Q1D, Q1D, Q1D);
60 const auto J = Reshape(j.Read(), Q1D, Q1D, Q1D, 3, 3, NE);
61 auto y = Reshape(op.Write(), Q1D, Q1D, Q1D, 3, 3, NE);
62
63 mfem::forall_3D(NE, Q1D, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
64 {
65 MFEM_FOREACH_THREAD_DIRECT(qx, x, Q1D)
66 {
67 MFEM_FOREACH_THREAD_DIRECT(qy, y, Q1D)
68 {
69 MFEM_FOREACH_THREAD_DIRECT(qz, z, Q1D)
70 {
71 const real_t J11 = J(qx, qy, qz, 0, 0, e);
72 const real_t J21 = J(qx, qy, qz, 1, 0, e);
73 const real_t J31 = J(qx, qy, qz, 2, 0, e);
74 const real_t J12 = J(qx, qy, qz, 0, 1, e);
75 const real_t J22 = J(qx, qy, qz, 1, 1, e);
76 const real_t J32 = J(qx, qy, qz, 2, 1, e);
77 const real_t J13 = J(qx, qy, qz, 0, 2, e);
78 const real_t J23 = J(qx, qy, qz, 1, 2, e);
79 const real_t J33 = J(qx, qy, qz, 2, 2, e);
80 const real_t cw = W(qx, qy, qz) * COEFF;
81 // adj(J)
82 const real_t A11 = (J22 * J33) - (J23 * J32);
83 const real_t A12 = (J32 * J13) - (J12 * J33);
84 const real_t A13 = (J12 * J23) - (J22 * J13);
85 const real_t A21 = (J31 * J23) - (J21 * J33);
86 const real_t A22 = (J11 * J33) - (J13 * J31);
87 const real_t A23 = (J21 * J13) - (J11 * J23);
88 const real_t A31 = (J21 * J32) - (J31 * J22);
89 const real_t A32 = (J31 * J12) - (J11 * J32);
90 const real_t A33 = (J11 * J22) - (J12 * J21);
91 // Store wq * Q * adj(J)
92 y(qx, qy, qz, 0, 0, e) = cw * A11;
93 y(qx, qy, qz, 0, 1, e) = cw * A12;
94 y(qx, qy, qz, 0, 2, e) = cw * A13;
95 y(qx, qy, qz, 1, 0, e) = cw * A21;
96 y(qx, qy, qz, 1, 1, e) = cw * A22;
97 y(qx, qy, qz, 1, 2, e) = cw * A23;
98 y(qx, qy, qz, 2, 0, e) = cw * A31;
99 y(qx, qy, qz, 2, 1, e) = cw * A32;
100 y(qx, qy, qz, 2, 2, e) = cw * A33;
101 }
102 }
103 }
104 });
105}
106
107static void PADivergenceSetup(const int dim,
108 const int Q1D,
109 const int NE,
110 const Array<real_t> &W,
111 const Vector &J,
112 const real_t COEFF,
113 Vector &op)
114{
115 if (dim == 1) { MFEM_ABORT("dim==1 not supported in PADivergenceSetup"); }
116 else if (dim == 2)
117 {
118 PADivergenceSetup2D(Q1D, NE, W, J, COEFF, op);
119 }
120 else if (dim == 3)
121 {
122 PADivergenceSetup3D(Q1D, NE, W, J, COEFF, op);
123 }
124}
125
126void VectorDivergenceIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
127 const FiniteElementSpace &test_fes)
128{
129 // Assumes tensor-product elements ordered by nodes
130 MFEM_ASSERT(trial_fes.GetOrdering() == Ordering::byNODES,
131 "PA Only supports Ordering::byNODES!");
132 auto *mesh = trial_fes.GetMesh();
133 const auto &trial_fe = *trial_fes.GetTypicalFE();
134 const auto &test_fe = *test_fes.GetTypicalFE();
135 const auto *ir = IntRule ? IntRule :
136 &GetRule(trial_fe, test_fe,
137 *mesh->GetTypicalElementTransformation());
138 const int dims = trial_fe.GetDim();
139 nq = ir->GetNPoints();
140 dim = mesh->Dimension();
141 ne = trial_fes.GetNE();
142 geom = mesh->GetGeometricFactors(*ir, GeometricFactors::JACOBIANS);
143 trial_maps = &trial_fe.GetDofToQuad(*ir, DofToQuad::TENSOR);
144 trial_dofs1D = trial_maps->ndof;
145 quad1D = trial_maps->nqpt;
146 test_maps = &test_fe.GetDofToQuad(*ir, DofToQuad::TENSOR);
147 test_dofs1D = test_maps->ndof;
148 MFEM_ASSERT(quad1D == test_maps->nqpt,
149 "PA requires test and trial space to have same number of "
150 "quadrature points!");
151 pa_data.SetSize(nq * dims * dims * ne, Device::GetMemoryType());
152
153 real_t coeff = 1.0;
154 if (Q)
155 {
156 auto *cQ = dynamic_cast<ConstantCoefficient *>(Q);
157 MFEM_VERIFY(cQ, "only ConstantCoefficient is supported!");
158 coeff = cQ->constant;
159 }
160 PADivergenceSetup(dim, quad1D, ne, ir->GetWeights(), geom->J, coeff, pa_data);
161}
162
163// PA Divergence Apply kernel
164void VectorDivergenceIntegrator::AddMultPA(const Vector &x, Vector &y) const
165{
166 using Kernels = VectorDivergenceAddMultPA;
167
168 static const auto specializations =
169 ( // 2D
170 Kernels::Specialization<2, 2, 2, 2>::Add(),
171 Kernels::Specialization<2, 2, 2, 3>::Add(),
172 Kernels::Specialization<2, 2, 2, 4>::Add(),
173 Kernels::Specialization<2, 3, 2, 3>::Add(),
174 Kernels::Specialization<2, 3, 2, 4>::Add(),
175 Kernels::Specialization<2, 3, 3, 3>::Add(),
176 Kernels::Specialization<2, 3, 3, 4>::Add(),
177 Kernels::Specialization<2, 3, 3, 5>::Add(),
178 Kernels::Specialization<2, 4, 3, 4>::Add(),
179 Kernels::Specialization<2, 4, 3, 5>::Add(),
180 Kernels::Specialization<2, 4, 4, 4>::Add(),
181 Kernels::Specialization<2, 4, 4, 5>::Add(),
182 Kernels::Specialization<2, 4, 4, 6>::Add(),
183 Kernels::Specialization<2, 5, 4, 5>::Add(),
184 Kernels::Specialization<2, 5, 4, 6>::Add(),
185 Kernels::Specialization<2, 5, 5, 5>::Add(),
186 Kernels::Specialization<2, 5, 5, 6>::Add(),
187 Kernels::Specialization<2, 5, 5, 7>::Add(),
188 // 3D
189 Kernels::Specialization<3, 2, 2, 3>::Add(),
190 Kernels::Specialization<3, 2, 2, 4>::Add(),
191 Kernels::Specialization<3, 2, 2, 6>::Add(),
192 Kernels::Specialization<3, 3, 2, 5>::Add(),
193 Kernels::Specialization<3, 3, 3, 4>::Add(),
194 Kernels::Specialization<3, 3, 3, 5>::Add(),
195 Kernels::Specialization<3, 3, 3, 7>::Add(),
196 Kernels::Specialization<3, 4, 4, 5>::Add(),
197 Kernels::Specialization<3, 4, 4, 6>::Add(),
198 Kernels::Specialization<3, 4, 4, 8>::Add(),
199 Kernels::Specialization<3, 5, 5, 6>::Add(),
200 Kernels::Specialization<3, 5, 5, 7>::Add(),
201 Kernels::Specialization<3, 5, 5, 9>::Add(),
202 true);
203 MFEM_CONTRACT_VAR(specializations);
204
205 Kernels::Run(dim, trial_dofs1D, test_dofs1D, quad1D, ne,
206 trial_maps->B, trial_maps->G, test_maps->Bt,
207 pa_data, x, y,
208 trial_dofs1D, test_dofs1D, quad1D);
209}
210
211// PA Divergence Apply kernel transpose
212void VectorDivergenceIntegrator::AddMultTransposePA(const Vector &x,
213 Vector &y) const
214{
215 using Kernels = VectorDivergenceAddMultTransposePA;
216
217 static const auto specializations =
218 ( // 2D
219 Kernels::Specialization<2, 2, 2, 2>::Add(),
220 Kernels::Specialization<2, 2, 2, 3>::Add(),
221 Kernels::Specialization<2, 2, 2, 4>::Add(),
222 Kernels::Specialization<2, 3, 2, 3>::Add(),
223 Kernels::Specialization<2, 3, 2, 4>::Add(),
224 Kernels::Specialization<2, 3, 3, 3>::Add(),
225 Kernels::Specialization<2, 3, 3, 4>::Add(),
226 Kernels::Specialization<2, 3, 3, 5>::Add(),
227 Kernels::Specialization<2, 4, 3, 4>::Add(),
228 Kernels::Specialization<2, 4, 3, 5>::Add(),
229 Kernels::Specialization<2, 4, 4, 4>::Add(),
230 Kernels::Specialization<2, 4, 4, 5>::Add(),
231 Kernels::Specialization<2, 4, 4, 6>::Add(),
232 Kernels::Specialization<2, 5, 4, 5>::Add(),
233 Kernels::Specialization<2, 5, 4, 6>::Add(),
234 Kernels::Specialization<2, 5, 5, 5>::Add(),
235 Kernels::Specialization<2, 5, 5, 6>::Add(),
236 Kernels::Specialization<2, 5, 5, 7>::Add(),
237 // 3D
238 Kernels::Specialization<3, 2, 2, 3>::Add(),
239 Kernels::Specialization<3, 2, 2, 4>::Add(),
240 Kernels::Specialization<3, 2, 2, 6>::Add(),
241 Kernels::Specialization<3, 3, 2, 5>::Add(),
242 Kernels::Specialization<3, 3, 3, 4>::Add(),
243 Kernels::Specialization<3, 3, 3, 5>::Add(),
244 Kernels::Specialization<3, 3, 3, 7>::Add(),
245 Kernels::Specialization<3, 4, 4, 5>::Add(),
246 Kernels::Specialization<3, 4, 4, 6>::Add(),
247 Kernels::Specialization<3, 4, 4, 8>::Add(),
248 Kernels::Specialization<3, 5, 5, 6>::Add(),
249 Kernels::Specialization<3, 5, 5, 7>::Add(),
250 Kernels::Specialization<3, 5, 5, 9>::Add(),
251 true);
252 MFEM_CONTRACT_VAR(specializations);
253
254 Kernels::Run(dim, trial_dofs1D, test_dofs1D, quad1D, ne,
255 trial_maps->Bt, trial_maps->Gt, test_maps->B,
256 pa_data, x, y,
257 trial_dofs1D, test_dofs1D, quad1D);
258}
259
260} // namespace mfem
bilininteg_vecdiv_pa.hpp
mfem::ConstantCoefficient
A coefficient that is constant across space and time.
Definition coefficient.hpp:88
mfem::Device::GetMemoryType
static MemoryType GetMemoryType()
(DEPRECATED) Equivalent to GetDeviceMemoryType().
Definition device.hpp:302
mfem::DofToQuad::G
Array< real_t > G
Gradients/divergences/curls of basis functions evaluated at quadrature points.
Definition fe_base.hpp:222
mfem::DofToQuad::TENSOR
@ TENSOR
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition fe_base.hpp:165
mfem::DofToQuad::B
Array< real_t > B
Basis functions evaluated at quadrature points.
Definition fe_base.hpp:201
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_base.hpp:186
mfem::DofToQuad::Gt
Array< real_t > Gt
Transpose of G.
Definition fe_base.hpp:229
mfem::DofToQuad::nqpt
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Definition fe_base.hpp:190
mfem::DofToQuad::Bt
Array< real_t > Bt
Transpose of B.
Definition fe_base.hpp:207
mfem::FiniteElementSpace
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:210
mfem::FiniteElementSpace::GetOrdering
Ordering::Type GetOrdering() const
Return the ordering method.
Definition fespace.hpp:852
mfem::FiniteElementSpace::GetNE
int GetNE() const
Returns number of elements in the mesh.
Definition fespace.hpp:867
mfem::FiniteElementSpace::GetMesh
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:639
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:3896
mfem::GeometricFactors::J
Vector J
Jacobians of the element transformations at all quadrature points.
Definition mesh.hpp:3158
mfem::Integrator::IntRule
const IntegrationRule * IntRule
Definition integrator.hpp:70
mfem::VectorDivergenceIntegrator::AssemblePA
void AssemblePA(const FiniteElementSpace &trial_fes, const FiniteElementSpace &test_fes) override
Definition bilininteg_vecdiv_pa.cpp:126
mfem::VectorDivergenceIntegrator::AddMultTransposePA
void AddMultTransposePA(const Vector &x, Vector &y) const override
Method for partially assembled transposed action.
Definition bilininteg_vecdiv_pa.cpp:212
mfem::VectorDivergenceIntegrator::GetRule
static const IntegrationRule & GetRule(const FiniteElement &trial_fe, const FiniteElement &test_fe, const ElementTransformation &Trans)
Definition bilininteg.cpp:2942
mfem::VectorDivergenceIntegrator::AddMultPA
void AddMultPA(const Vector &x, Vector &y) const override
Method for partially assembled action.
Definition bilininteg_vecdiv_pa.cpp:164
mfem::Vector
Vector data type.
Definition vector.hpp:82
mfem::Vector::SetSize
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:633
dim
int dim
Definition ex24.cpp:53
real_t
mfem::real_t real_t
Definition mtop_solvers.hpp:19
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:138
mfem::forall_2D
void forall_2D(int N, int X, int Y, lambda &&body)
Definition forall.hpp:1220
mfem::forall_3D
void forall_3D(int N, int X, int Y, int Z, lambda &&body)
Definition forall.hpp:1244
mfem::real_t
float real_t
Definition config.hpp:46
mfem::SchrodingerBaseKernels
Base class for Schrodinger solver kernels.
Definition schrodinger_flow.hpp:242
Generated on Tue Sep 1 2026 16:39:57 for MFEM by doxygen 1.11.0