MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
coefficient.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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#ifndef MFEM_LIBCEED_COEFF
13#define MFEM_LIBCEED_COEFF
14
15#ifdef MFEM_USE_CEED
16
22#include "../../gridfunc.hpp"
23#include "../../qfunction.hpp"
24#include "util.hpp"
25#include "ceed.hpp"
26
27namespace mfem
28{
29
30class Mesh;
31class IntegrationRule;
32class Coefficient;
33class VectorCoefficient;
34class GridFunction;
35
36namespace ceed
37{
38
40{
41 const int ncomp;
42 Coefficient(int ncomp_) : ncomp(ncomp_) { }
43 virtual bool IsConstant() const { return true; }
44 virtual ~Coefficient() { }
45};
46
48{
49 CeedVector coeffVector = nullptr;
50 const CeedEvalMode emode;
51 VariableCoefficient(int ncomp_, CeedEvalMode emode_)
52 : Coefficient(ncomp_), emode(emode_) { }
53 virtual bool IsConstant() const override { return false; }
55 {
56 CeedVectorDestroy(&coeffVector);
57 }
58};
59
61{
63 CeedBasis basis;
64 CeedElemRestriction restr;
66 : VariableCoefficient(gf_.VectorDim(), CEED_EVAL_INTERP),
67 gf(gf_)
68 {
70 }
71};
72
74{
76 CeedElemRestriction restr;
77 QuadCoefficient(int ncomp_) : VariableCoefficient(ncomp_, CEED_EVAL_NONE) { }
78};
79
80/** @brief Initializes an mfem::ceed::Coefficient @a coeff_ptr from an
81 mfem::Coefficient @a Q, an mfem::Mesh @a mesh, and an mfem::IntegrationRule
82 @a ir.
83
84 @param[in] Q is the coefficient from the `Integrator`.
85 @param[in] mesh is the mesh.
86 @param[in] ir is the integration rule.
87 @param[out] coeff_ptr is the structure to store the coefficient for the
88 `CeedOperator`.
89 @param[out] ctx is the Context associated to the QFunction. */
90template <typename Context>
92 const mfem::IntegrationRule &ir,
93 Coefficient*& coeff_ptr, Context &ctx)
94{
95 if ( Q == nullptr )
96 {
97 Coefficient *ceedCoeff = new Coefficient(1);
98 ctx.coeff = 1.0;
99 coeff_ptr = ceedCoeff;
100 }
101 else if (ConstantCoefficient *const_coeff =
102 dynamic_cast<ConstantCoefficient*>(Q))
103 {
104 Coefficient *ceedCoeff = new Coefficient(1);
105 ctx.coeff = const_coeff->constant;
106 coeff_ptr = ceedCoeff;
107 }
108 else if (GridFunctionCoefficient* gf_coeff =
109 dynamic_cast<GridFunctionCoefficient*>(Q))
110 {
111 GridCoefficient *ceedCoeff =
112 new GridCoefficient(*gf_coeff->GetGridFunction());
113 coeff_ptr = ceedCoeff;
114 }
115 else if (QuadratureFunctionCoefficient *cQ =
116 dynamic_cast<QuadratureFunctionCoefficient*>(Q))
117 {
118 QuadCoefficient *ceedCoeff = new QuadCoefficient(1);
119 const int ne = mesh.GetNE();
120 const int nq = ir.GetNPoints();
121 const mfem::QuadratureFunction &qFun = cQ->GetQuadFunction();
122 MFEM_VERIFY(qFun.Size() == nq * ne,
123 "Incompatible QuadratureFunction dimension \n");
124
125 MFEM_VERIFY(&ir == &qFun.GetSpace()->GetIntRule(0),
126 "IntegrationRule used within integrator and in"
127 " QuadratureFunction appear to be different");
128 qFun.Read();
129 ceedCoeff->coeff.MakeRef(const_cast<mfem::QuadratureFunction &>(qFun),0);
130 InitVector(ceedCoeff->coeff, ceedCoeff->coeffVector);
131 coeff_ptr = ceedCoeff;
132 }
133 else
134 {
135 QuadCoefficient *ceedCoeff = new QuadCoefficient(1);
136 const int ne = mesh.GetNE();
137 const int nq = ir.GetNPoints();
138 ceedCoeff->coeff.SetSize(nq * ne);
139 auto C = Reshape(ceedCoeff->coeff.HostWrite(), nq, ne);
140 for (int e = 0; e < ne; ++e)
141 {
143 for (int q = 0; q < nq; ++q)
144 {
145 C(q,e) = Q->Eval(T, ir.IntPoint(q));
146 }
147 }
148 InitVector(ceedCoeff->coeff, ceedCoeff->coeffVector);
149 coeff_ptr = ceedCoeff;
150 }
151}
152
153
154/** @brief Initializes an mfem::ceed::Coefficient @a coeff_ptr from an
155 mfem::VectorCoefficient @a VQ, an mfem::Mesh @a mesh, and an
156 mfem::IntegrationRule @a ir.
157
158 @param[in] VQ is the vector coefficient from the `Integrator`.
159 @param[in] mesh is the mesh.
160 @param[in] ir is the integration rule.
161 @param[out] coeff_ptr is the structure to store the coefficient for the
162 `CeedOperator`.
163 @param[out] ctx is the Context associated to the QFunction. */
164template <typename Context>
166 const mfem::IntegrationRule &ir,
167 Coefficient *&coeff_ptr, Context &ctx)
168{
169 if (VectorConstantCoefficient *const_coeff =
170 dynamic_cast<VectorConstantCoefficient*>(VQ))
171 {
172 const int vdim = const_coeff->GetVDim();
173 const mfem::Vector &val = const_coeff->GetVec();
174 Coefficient *ceedCoeff = new Coefficient(vdim);
175 for (int i = 0; i < vdim; i++)
176 {
177 ctx.coeff[i] = val[i];
178 }
179 coeff_ptr = ceedCoeff;
180 }
181 else if (VectorGridFunctionCoefficient* vgf_coeff =
182 dynamic_cast<VectorGridFunctionCoefficient*>(VQ))
183 {
184 GridCoefficient *ceedCoeff =
185 new GridCoefficient(*vgf_coeff->GetGridFunction());
186 coeff_ptr = ceedCoeff;
187 }
189 dynamic_cast<VectorQuadratureFunctionCoefficient*>(VQ))
190 {
191 QuadCoefficient *ceedCoeff = new QuadCoefficient(cQ->GetVDim());
192 const int dim = mesh.Dimension();
193 const int ne = mesh.GetNE();
194 const int nq = ir.GetNPoints();
195 const mfem::QuadratureFunction &qFun = cQ->GetQuadFunction();
196 MFEM_VERIFY(qFun.Size() == dim * nq * ne,
197 "Incompatible QuadratureFunction dimension \n");
198
199 MFEM_VERIFY(&ir == &qFun.GetSpace()->GetIntRule(0),
200 "IntegrationRule used within integrator and in"
201 " QuadratureFunction appear to be different");
202 qFun.Read();
203 ceedCoeff->coeff.MakeRef(const_cast<mfem::QuadratureFunction &>(qFun),0);
204 InitVector(ceedCoeff->coeff, ceedCoeff->coeffVector);
205 coeff_ptr = ceedCoeff;
206 }
207 else
208 {
209 const int dim = mesh.Dimension();
210 QuadCoefficient *ceedCoeff = new QuadCoefficient(dim);
211 const int ne = mesh.GetNE();
212 const int nq = ir.GetNPoints();
213 ceedCoeff->coeff.SetSize(dim * nq * ne);
214 auto C = Reshape(ceedCoeff->coeff.HostWrite(), dim, nq, ne);
216 for (int e = 0; e < ne; ++e)
217 {
219 VQ->Eval(Q_ir, T, ir);
220 for (int q = 0; q < nq; ++q)
221 {
222 for (int i = 0; i < dim; ++i)
223 {
224 C(i,q,e) = Q_ir(i,q);
225 }
226 }
227 }
228 InitVector(ceedCoeff->coeff, ceedCoeff->coeffVector);
229 coeff_ptr = ceedCoeff;
230 }
231}
232
233/** @brief Initializes an mfem::ceed::Coefficient @a coeff_ptr from an
234 mfem::Coefficient @a Q, an mfem::Mesh @a mesh, and an mfem::IntegrationRule
235 @a ir for the elements given by the indices @a indices.
236
237 @param[in] Q is the coefficient from the `Integrator`.
238 @param[in] mesh is the mesh.
239 @param[in] ir is the integration rule.
240 @param[in] nelem The number of elements.
241 @param[in] indices The indices of the elements of same type in the
242 `FiniteElementSpace`.
243 @param[out] coeff_ptr is the structure to store the coefficient for the
244 `CeedOperator`.
245 @param[out] ctx is the Context associated to the QFunction. */
246template <typename Context>
248 const mfem::IntegrationRule &ir,
249 int nelem,
250 const int* indices,
251 Coefficient*& coeff_ptr, Context &ctx)
252{
253 if ( Q == nullptr )
254 {
255 Coefficient *ceedCoeff = new Coefficient(1);
256 ctx.coeff = 1.0;
257 coeff_ptr = ceedCoeff;
258 }
259 else if (ConstantCoefficient *const_coeff =
260 dynamic_cast<ConstantCoefficient*>(Q))
261 {
262 Coefficient *ceedCoeff = new Coefficient(1);
263 ctx.coeff = const_coeff->constant;
264 coeff_ptr = ceedCoeff;
265 }
266 else if (GridFunctionCoefficient* gf_coeff =
267 dynamic_cast<GridFunctionCoefficient*>(Q))
268 {
269 GridCoefficient *ceedCoeff =
270 new GridCoefficient(*gf_coeff->GetGridFunction());
271 coeff_ptr = ceedCoeff;
272 }
273 else if (QuadratureFunctionCoefficient *cQ =
274 dynamic_cast<QuadratureFunctionCoefficient*>(Q))
275 {
276 QuadCoefficient *ceedCoeff = new QuadCoefficient(1);
277 const int ne = mesh.GetNE();
278 const int nq = ir.GetNPoints();
279 const mfem::QuadratureFunction &qFun = cQ->GetQuadFunction();
280 MFEM_VERIFY(qFun.Size() == nq * ne,
281 "Incompatible QuadratureFunction dimension \n");
282
283 MFEM_VERIFY(&ir == &qFun.GetSpace()->GetIntRule(0),
284 "IntegrationRule used within integrator and in"
285 " QuadratureFunction appear to be different");
286 ceedCoeff->coeff.SetSize(nq * nelem);
287 Memory<int> m_indices((int*)indices, nelem, false);
288 auto in = Reshape(qFun.Read(), nq, ne);
289 auto d_indices = Read(m_indices, nelem);
290 auto out = Reshape(ceedCoeff->coeff.Write(), nq, nelem);
291 mfem::forall(nelem * nq, [=] MFEM_HOST_DEVICE (int i)
292 {
293 const int q = i%nq;
294 const int sub_e = i/nq;
295 const int e = d_indices[sub_e];
296 out(q, sub_e) = in(q, e);
297 });
298 m_indices.DeleteDevice();
299 InitVector(ceedCoeff->coeff, ceedCoeff->coeffVector);
300 coeff_ptr = ceedCoeff;
301 }
302 else
303 {
304 QuadCoefficient *ceedCoeff = new QuadCoefficient(1);
305 const int nq = ir.GetNPoints();
306 ceedCoeff->coeff.SetSize(nq * nelem);
307 auto C = Reshape(ceedCoeff->coeff.HostWrite(), nq, nelem);
308 for (int i = 0; i < nelem; ++i)
309 {
310 const int e = indices[i];
312 for (int q = 0; q < nq; ++q)
313 {
314 C(q, i) = Q->Eval(T, ir.IntPoint(q));
315 }
316 }
317 InitVector(ceedCoeff->coeff, ceedCoeff->coeffVector);
318 coeff_ptr = ceedCoeff;
319 }
320}
321
322
323/** @brief Initializes an mfem::ceed::Coefficient @a coeff_ptr from an
324 mfem::VectorCoefficient @a Q, an mfem::Mesh @a mesh, and an
325 mfem::IntegrationRule @a ir for the elements given by the indices @a indices.
326
327 @param[in] VQ is the vector coefficient from the `Integrator`.
328 @param[in] mesh is the mesh.
329 @param[in] ir is the integration rule.
330 @param[in] nelem The number of elements.
331 @param[in] indices The indices of the elements of same type in the
332 `FiniteElementSpace`.
333 @param[out] coeff_ptr is the structure to store the coefficient for the
334 `CeedOperator`.
335 @param[out] ctx is the Context associated to the QFunction. */
336template <typename Context>
338 const mfem::IntegrationRule &ir,
339 int nelem,
340 const int* indices,
341 Coefficient *&coeff_ptr, Context &ctx)
342{
343 if (VectorConstantCoefficient *const_coeff =
344 dynamic_cast<VectorConstantCoefficient*>(VQ))
345 {
346 const int vdim = const_coeff->GetVDim();
347 const mfem::Vector &val = const_coeff->GetVec();
348 Coefficient *ceedCoeff = new Coefficient(vdim);
349 for (int i = 0; i < vdim; i++)
350 {
351 ctx.coeff[i] = val[i];
352 }
353 coeff_ptr = ceedCoeff;
354 }
355 else if (VectorGridFunctionCoefficient* vgf_coeff =
356 dynamic_cast<VectorGridFunctionCoefficient*>(VQ))
357 {
358 GridCoefficient *ceedCoeff =
359 new GridCoefficient(*vgf_coeff->GetGridFunction());
360 coeff_ptr = ceedCoeff;
361 }
363 dynamic_cast<VectorQuadratureFunctionCoefficient*>(VQ))
364 {
365 QuadCoefficient *ceedCoeff = new QuadCoefficient(cQ->GetVDim());
366 const int dim = mesh.Dimension();
367 const int ne = mesh.GetNE();
368 const int nq = ir.GetNPoints();
369 const mfem::QuadratureFunction &qFun = cQ->GetQuadFunction();
370 MFEM_VERIFY(qFun.Size() == dim * nq * ne,
371 "Incompatible QuadratureFunction dimension \n");
372
373 MFEM_VERIFY(&ir == &qFun.GetSpace()->GetIntRule(0),
374 "IntegrationRule used within integrator and in"
375 " QuadratureFunction appear to be different");
376 ceedCoeff->coeff.SetSize(dim * nq * nelem);
377 Memory<int> m_indices((int*)indices, nelem, false);
378 auto in = Reshape(qFun.Read(), dim, nq, ne);
379 auto d_indices = Read(m_indices, nelem);
380 auto out = Reshape(ceedCoeff->coeff.Write(), dim, nq, nelem);
381 mfem::forall(nelem * nq, [=] MFEM_HOST_DEVICE (int i)
382 {
383 const int q = i%nq;
384 const int sub_e = i/nq;
385 const int e = d_indices[sub_e];
386 for (int d = 0; d < dim; d++)
387 {
388 out(d, q, sub_e) = in(d, q, e);
389 }
390 });
391 m_indices.DeleteDevice();
392 InitVector(ceedCoeff->coeff, ceedCoeff->coeffVector);
393 coeff_ptr = ceedCoeff;
394 }
395 else
396 {
397 const int dim = mesh.Dimension();
398 QuadCoefficient *ceedCoeff = new QuadCoefficient(dim);
399 const int nq = ir.GetNPoints();
400 ceedCoeff->coeff.SetSize(dim * nq * nelem);
401 auto C = Reshape(ceedCoeff->coeff.HostWrite(), dim, nq, nelem);
403 for (int i = 0; i < nelem; ++i)
404 {
405 const int e = indices[i];
407 VQ->Eval(Q_ir, T, ir);
408 for (int q = 0; q < nq; ++q)
409 {
410 for (int d = 0; d < dim; ++d)
411 {
412 C(d, q, i) = Q_ir(d, q);
413 }
414 }
415 }
416 InitVector(ceedCoeff->coeff, ceedCoeff->coeffVector);
417 coeff_ptr = ceedCoeff;
418 }
419}
420
421template <typename Coeff, typename Context>
423 const mfem::IntegrationRule &ir, int nelem,
424 const int* indices, Coefficient *&coeff_ptr, Context &ctx)
425{
426 if (indices)
427 {
428 InitCoefficientWithIndices(Q, mesh, ir, nelem, indices, coeff_ptr, ctx);
429 }
430 else
431 {
432 InitCoefficient(Q, mesh, ir, coeff_ptr, ctx);
433 }
434}
435
436} // namespace ceed
437
438} // namespace mfem
439
440#endif
441
442#endif // MFEM_LIBCEED_COEFF
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
A coefficient that is constant across space and time.
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition intrules.hpp:256
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition intrules.hpp:259
Class used by MFEM to store pointers to host and/or device memory.
void DeleteDevice(bool copy_to_host=true)
Delete the device pointer, if owned. If copy_to_host is true and the data is valid only on device,...
Mesh data type.
Definition mesh.hpp:56
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void GetElementTransformation(int i, IsoparametricTransformation *ElTr) const
Builds the transformation defining the i-th element in ElTr. ElTr must be allocated in advance and wi...
Definition mesh.cpp:357
Quadrature function coefficient which requires that the quadrature rules used for this coefficient be...
Represents values or vectors of values at quadrature points on a mesh.
Definition qfunction.hpp:24
QuadratureSpaceBase * GetSpace()
Get the associated QuadratureSpaceBase object.
Definition qfunction.hpp:82
const IntegrationRule & GetIntRule(int idx) const
Return the IntegrationRule associated with entity idx.
Definition qspace.hpp:81
Base class for vector Coefficients that optionally depend on time and space.
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
Vector coefficient that is constant in space and time.
Vector coefficient defined by a vector GridFunction.
Vector quadrature function coefficient which requires that the quadrature rules used for this vector ...
Vector data type.
Definition vector.hpp:80
virtual const real_t * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:474
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
virtual real_t * HostWrite()
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
Definition vector.hpp:486
virtual real_t * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:482
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition vector.hpp:602
int dim
Definition ex24.cpp:53
void InitVector(const mfem::Vector &v, CeedVector &cv)
Initialize a CeedVector from an mfem::Vector.
Definition util.cpp:75
void InitCoefficient(mfem::Coefficient *Q, mfem::Mesh &mesh, const mfem::IntegrationRule &ir, Coefficient *&coeff_ptr, Context &ctx)
Initializes an mfem::ceed::Coefficient coeff_ptr from an mfem::Coefficient Q, an mfem::Mesh mesh,...
void InitCoefficientWithIndices(mfem::Coefficient *Q, mfem::Mesh &mesh, const mfem::IntegrationRule &ir, int nelem, const int *indices, Coefficient *&coeff_ptr, Context &ctx)
Initializes an mfem::ceed::Coefficient coeff_ptr from an mfem::Coefficient Q, an mfem::Mesh mesh,...
const T * Read(const Memory< T > &mem, int size, bool on_dev=true)
Get a pointer for read access to mem with the mfem::Device's DeviceMemoryClass, if on_dev = true,...
Definition device.hpp:320
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition globals.hpp:66
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
void forall(int N, lambda &&body)
Definition forall.hpp:754
struct s_NavierContext ctx
virtual bool IsConstant() const
CeedElemRestriction restr
GridCoefficient(const mfem::GridFunction &gf_)
const mfem::GridFunction & gf
CeedElemRestriction restr
virtual bool IsConstant() const override
VariableCoefficient(int ncomp_, CeedEvalMode emode_)