MFEM  v4.6.0 Finite element discretization library
coefficient.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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
17 #include "../../../general/forall.hpp"
18 #include "../../../config/config.hpp"
19 #include "../../../linalg/vector.hpp"
20 #include "../../../linalg/dtensor.hpp"
21 #include "../../../mesh/mesh.hpp"
22 #include "../../gridfunc.hpp"
23 #include "../../qfunction.hpp"
24 #include "util.hpp"
25 #include "ceed.hpp"
26
27 namespace mfem
28 {
29
30 class Mesh;
31 class IntegrationRule;
32 class Coefficient;
33 class VectorCoefficient;
34 class GridFunction;
35
36 namespace 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. */
90 template <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. */
164 template <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);
215  mfem::DenseMatrix Q_ir;
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. */
246 template <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. */
336 template <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);
402  mfem::DenseMatrix Q_ir;
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
421 template <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
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:253
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:96
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
Base class for vector Coefficients that optionally depend on time and space.
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
void InitVector(const mfem::Vector &v, CeedVector &cv)
Initialize a CeedVector from an mfem::Vector.
Definition: util.cpp:75
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 ...
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
virtual bool IsConstant() const
Definition: coefficient.hpp:43
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
CeedElemRestriction restr
Definition: coefficient.hpp:76
Vector quadrature function coefficient which requires that the quadrature rules used for this vector ...
Vector coefficient that is constant in space and time.
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
virtual double * HostWrite()
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:465
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:453
const mfem::GridFunction & gf
Definition: coefficient.hpp:62
GridCoefficient(const mfem::GridFunction &gf_)
Definition: coefficient.hpp:65
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:256
Quadrature function coefficient which requires that the quadrature rules used for this coefficient be...
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...
Definition: coefficient.hpp:91
virtual double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:461
VariableCoefficient(int ncomp_, CeedEvalMode emode_)
Definition: coefficient.hpp:51
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&#39;s DeviceMemoryClass, if on_dev = true...
Definition: device.hpp:319
void forall(int N, lambda &&body)
Definition: forall.hpp:742
CeedElemRestriction restr
Definition: coefficient.hpp:64
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
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
const IntegrationRule & GetIntRule(int idx) const
Return the IntegrationRule associated with entity idx.
Definition: qspace.hpp:73
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...
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
QuadratureSpaceBase * GetSpace()
Get the associated QuadratureSpaceBase object.
Definition: qfunction.hpp:82
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:355
virtual bool IsConstant() const override
Definition: coefficient.hpp:53
int dim
Definition: ex24.cpp:53
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
Vector data type.
Definition: vector.hpp:58
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...
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition: vector.hpp:581
Vector coefficient defined by a vector GridFunction.
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
Represents values or vectors of values at quadrature points on a mesh.
Definition: qfunction.hpp:23