40 MFEM_VERIFY(vdim == fes.
GetVDim(),
"vdim != fes.GetVDim()");
41 MFEM_VERIFY(vdim == mesh->
Dimension(),
"vdim != dim");
49 const int nq = ir->GetNPoints();
57 if (!(
dim == 2 ||
dim == 3)) { MFEM_ABORT(
"Dimension not supported."); }
69 MFEM_VERIFY(
VQ->
GetVDim() == vdim,
"VQ vdim vs. vdim error");
74 MFEM_VERIFY(
MQ->
GetVDim() == vdim,
"MQ dimension vs. vdim error");
75 MFEM_VERIFY(coeff.
Size() == (vdim*vdim) *
ne * nq,
"MQ size error");
82 const bool matrix_coeff =
coeff_vdim == vdim * vdim;
83 MFEM_VERIFY(const_coeff + vector_coeff + matrix_coeff == 1,
"");
87 const auto w_r = ir->GetWeights().Read();
91 const auto W =
Reshape(w_r, q1d, q1d);
98 MFEM_FOREACH_THREAD(qy, y, q1d)
100 MFEM_FOREACH_THREAD(qx, x, q1d)
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;
119 const auto W =
Reshape(w_r, q1d, q1d, q1d);
126 MFEM_FOREACH_THREAD(qz, z, q1d)
128 MFEM_FOREACH_THREAD(qy, y, q1d)
130 MFEM_FOREACH_THREAD(qx, x, q1d)
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;
163 MFEM_ABORT(
"Unknown VectorMassIntegrator::AssemblePA kernel for"
164 <<
" dim:" <<
dim <<
", vdim:" << vdim <<
", sdim:" << sdim);
174 static const auto vector_mass_kernel_specializations =
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(),
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(),
200 MFEM_CONTRACT_VAR(vector_mass_kernel_specializations);
208template <const
int T_D1D = 0, const
int T_Q1D = 0>
209static void PAVectorMassAssembleDiagonal2D(
const int NE,
212 const int d1d = 0,
const int q1d = 0)
214 constexpr int VDIM = 2;
215 const int D1D = T_D1D ? T_D1D : d1d;
216 const int Q1D = T_Q1D ? T_Q1D : q1d;
219 const auto B =
Reshape(
b.Read(), Q1D, D1D);
220 const auto D =
Reshape(pa_data.
Read(), Q1D, Q1D, NE);
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;
230 real_t temp[max_Q1D][max_D1D];
231 for (
int qx = 0; qx < Q1D; ++qx)
233 for (
int dy = 0; dy < D1D; ++dy)
236 for (
int qy = 0; qy < Q1D; ++qy)
238 temp[qx][dy] += B(qy, dy) * B(qy, dy) * D(qx, qy, e);
242 for (
int dy = 0; dy < D1D; ++dy)
244 for (
int dx = 0; dx < D1D; ++dx)
247 for (
int qx = 0; qx < Q1D; ++qx)
249 temp1 += B(qx, dx) * B(qx, dx) * temp[qx][dy];
251 Y(dx, dy, 0, e) = temp1;
252 Y(dx, dy, 1, e) = temp1;
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)
264 constexpr int VDIM = 3;
265 const int D1D = T_D1D ? T_D1D : d1d;
266 const int Q1D = T_Q1D ? T_Q1D : 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);
275 const int D1D = T_D1D ? T_D1D : d1d;
276 const int Q1D = T_Q1D ? T_Q1D : q1d;
278 constexpr int max_D1D = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
279 constexpr int max_Q1D = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
281 real_t temp[max_Q1D][max_Q1D][max_D1D];
282 for (
int qx = 0; qx < Q1D; ++qx)
284 for (
int qy = 0; qy < Q1D; ++qy)
286 for (
int dz = 0; dz < D1D; ++dz)
288 temp[qx][qy][dz] = 0.0;
289 for (
int qz = 0; qz < Q1D; ++qz)
292 B(qz, dz) * B(qz, dz) * D(qx, qy, qz, e);
297 real_t temp2[max_Q1D][max_D1D][max_D1D];
298 for (
int qx = 0; qx < Q1D; ++qx)
300 for (
int dz = 0; dz < D1D; ++dz)
302 for (
int dy = 0; dy < D1D; ++dy)
304 temp2[qx][dy][dz] = 0.0;
305 for (
int qy = 0; qy < Q1D; ++qy)
308 B(qy, dy) * B(qy, dy) * temp[qx][qy][dz];
313 for (
int dz = 0; dz < D1D; ++dz)
315 for (
int dy = 0; dy < D1D; ++dy)
317 for (
int dx = 0; dx < D1D; ++dx)
320 for (
int qx = 0; qx < Q1D; ++qx)
322 temp3 += B(qx, dx) * B(qx, dx) * temp2[qx][dy][dz];
324 Y(dx, dy, dz, 0, e) = temp3;
325 Y(dx, dy, dz, 1, e) = temp3;
326 Y(dx, dy, dz, 2, e) = temp3;
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,
341 return PAVectorMassAssembleDiagonal2D(NE, B, pa_data, diag, D1D, Q1D);
345 return PAVectorMassAssembleDiagonal3D(NE, B, pa_data, diag, D1D, Q1D);
347 MFEM_ABORT(
"Dimension not implemented.");
355 MFEM_VERIFY(
coeff_vdim == 1,
"coeff_vdim != 1");
356 MFEM_VERIFY(!
VQ && !
MQ,
"VQ and MQ not supported");
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...
@ TENSOR
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Array< real_t > B
Basis functions evaluated at quadrature points.
int ndof
Number of degrees of freedom = number of basis functions. When mode is TENSOR, this is the 1D number.
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
bool IsVariableOrder() const
Returns true if the space contains elements of varying polynomial orders.
Mesh * GetMesh() const
Returns the mesh.
int GetVDim() const
Returns the vector dimension of the finite element space.
const FiniteElement * GetTypicalFE() const
Return GetFE(0) if the local mesh is not empty; otherwise return a typical FE based on the Geometry t...
Abstract class for all finite elements.
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...
Vector J
Jacobians of the element transformations at all quadrature points.
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.
int GetNE() const
Returns number of elements.
int Dimension() const
Dimension of the reference space used within the elements.
ElementTransformation * GetTypicalElementTransformation()
If the local mesh is not empty return GetElementTransformation(0); otherwise, return the identity tra...
int SpaceDimension() const
Dimension of the physical space containing the mesh.
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.
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
Class representing the storage layout of a QuadratureFunction.
int GetVDim()
Returns dimension of the vector.
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.
void AssemblePA(const FiniteElementSpace &fes) override
Method defining partial assembly.
const GeometricFactors * geom
Not owned.
virtual const real_t * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
virtual real_t * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
int Size() const
Returns the size of the vector.
void SetSize(int s)
Resize the vector to size s.
virtual real_t * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
void GetDiagonal(mfem::Vector &diag) const
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).
Represent a MassIntegrator with AssemblyLevel::Partial using libCEED.
MFEM_HOST_DEVICE DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims... dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
bool DeviceCanUseCeed()
Function that determines if a CEED kernel should be used, based on the current mfem::Device configura...
void forall_2D(int N, int X, int Y, lambda &&body)
void forall_3D(int N, int X, int Y, int Z, lambda &&body)
MemoryType
Memory types supported by MFEM.
void forall(int N, lambda &&body)
static const DeviceDofQuadLimits & Get()
Return a const reference to the DeviceDofQuadLimits singleton.