MFEM  v4.5.2
Finite element discretization library
qfunction.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_QFUNCTION
13 #define MFEM_QFUNCTION
14 
15 #include "../config/config.hpp"
16 #include "qspace.hpp"
17 #include "gridfunc.hpp"
18 
19 namespace mfem
20 {
21 
22 /// Represents values or vectors of values at quadrature points on a mesh.
23 class QuadratureFunction : public Vector
24 {
25 protected:
26  QuadratureSpaceBase *qspace; ///< Associated QuadratureSpaceBase object.
27  bool own_qspace; ///< Does this own the associated QuadratureSpaceBase?
28  int vdim; ///< Vector dimension.
29 
30 public:
31  /// Default constructor, results in an empty vector.
32  QuadratureFunction() : qspace(nullptr), own_qspace(false), vdim(0) { }
33 
34  /// Create a QuadratureFunction based on the given QuadratureSpaceBase.
35  /** The QuadratureFunction does not assume ownership of the
36  QuadratureSpaceBase.
37  @note The Vector data is not initialized. */
38  QuadratureFunction(QuadratureSpaceBase &qspace_, int vdim_ = 1)
39  : Vector(vdim_*qspace_.GetSize()),
40  qspace(&qspace_), own_qspace(false), vdim(vdim_)
41  { }
42 
43  /// Create a QuadratureFunction based on the given QuadratureSpaceBase.
44  /** The QuadratureFunction does not assume ownership of the
45  QuadratureSpaceBase.
46  @warning @a qspace_ may not be NULL. */
47  QuadratureFunction(QuadratureSpaceBase *qspace_, int vdim_ = 1)
48  : QuadratureFunction(*qspace_, vdim_) { }
49 
50  /** @brief Create a QuadratureFunction based on the given QuadratureSpaceBase,
51  using the external (host) data, @a qf_data. */
52  /** The QuadratureFunction does not assume ownership of the
53  QuadratureSpaceBase or the external data.
54  @warning @a qspace_ may not be NULL.
55  @note @a qf_data must be a valid **host** pointer (see the constructor
56  Vector::Vector(double *, int)). */
57  QuadratureFunction(QuadratureSpaceBase *qspace_, double *qf_data, int vdim_ = 1)
58  : Vector(qf_data, vdim_*qspace_->GetSize()),
59  qspace(qspace_), own_qspace(false), vdim(vdim_) { }
60 
61  /** @brief Copy constructor. The QuadratureSpace ownership flag, #own_qspace,
62  in the new object is set to false. */
64  : QuadratureFunction(*orig.qspace, orig.vdim)
65  {
66  Vector::operator=(orig);
67  }
68 
69  /// Read a QuadratureFunction from the stream @a in.
70  /** The QuadratureFunction assumes ownership of the read QuadratureSpace. */
71  QuadratureFunction(Mesh *mesh, std::istream &in);
72 
73  /// Get the vector dimension.
74  int GetVDim() const { return vdim; }
75 
76  /// Set the vector dimension, updating the size by calling Vector::SetSize().
77  void SetVDim(int vdim_)
78  { vdim = vdim_; SetSize(vdim*qspace->GetSize()); }
79 
80  /// Get the associated QuadratureSpaceBase object.
82 
83  /// Get the associated QuadratureSpaceBase object (const version).
84  const QuadratureSpaceBase *GetSpace() const { return qspace; }
85 
86  /// Change the QuadratureSpaceBase and optionally the vector dimension.
87  /** If the new QuadratureSpaceBase is different from the current one, the
88  QuadratureFunction will not assume ownership of the new space; otherwise,
89  the ownership flag remains the same.
90 
91  If the new vector dimension @a vdim_ < 0, the vector dimension remains
92  the same.
93 
94  The data size is updated by calling Vector::SetSize(). */
95  inline void SetSpace(QuadratureSpaceBase *qspace_, int vdim_ = -1);
96 
97  /** @brief Change the QuadratureSpaceBase, the data array, and optionally the
98  vector dimension. */
99  /** If the new QuadratureSpaceBase is different from the current one, the
100  QuadratureFunction will not assume ownership of the new space; otherwise,
101  the ownership flag remains the same.
102 
103  If the new vector dimension @a vdim_ < 0, the vector dimension remains
104  the same.
105 
106  The data array is replaced by calling Vector::NewDataAndSize(). */
107  inline void SetSpace(QuadratureSpaceBase *qspace_, double *qf_data,
108  int vdim_ = -1);
109 
110  /// Get the QuadratureSpaceBase ownership flag.
111  bool OwnsSpace() { return own_qspace; }
112 
113  /// Set the QuadratureSpaceBase ownership flag.
114  void SetOwnsSpace(bool own) { own_qspace = own; }
115 
116  /// Set this equal to a constant value.
117  QuadratureFunction &operator=(double value);
118 
119  /// Copy the data from @a v.
120  /** The size of @a v must be equal to the size of the associated
121  QuadratureSpaceBase #qspace times the QuadratureFunction vector
122  dimension i.e. QuadratureFunction::Size(). */
124 
125  /// Evaluate a grid function at each quadrature point.
126  void ProjectGridFunction(const GridFunction &gf);
127 
128  /// Return all values associated with mesh element @a idx in a Vector.
129  /** The result is stored in the Vector @a values as a reference to the
130  global values.
131 
132  Inside the Vector @a values, the index `i+vdim*j` corresponds to the
133  `i`-th vector component at the `j`-th quadrature point.
134  */
135  inline void GetValues(int idx, Vector &values);
136 
137  /// Return all values associated with mesh element @a idx in a Vector.
138  /** The result is stored in the Vector @a values as a copy of the
139  global values.
140 
141  Inside the Vector @a values, the index `i+vdim*j` corresponds to the
142  `i`-th vector component at the `j`-th quadrature point.
143  */
144  inline void GetValues(int idx, Vector &values) const;
145 
146  /// Return the quadrature function values at an integration point.
147  /** The result is stored in the Vector @a values as a reference to the
148  global values. */
149  inline void GetValues(int idx, const int ip_num, Vector &values);
150 
151  /// Return the quadrature function values at an integration point.
152  /** The result is stored in the Vector @a values as a copy to the
153  global values. */
154  inline void GetValues(int idx, const int ip_num, Vector &values) const;
155 
156  /// Return all values associated with mesh element @a idx in a DenseMatrix.
157  /** The result is stored in the DenseMatrix @a values as a reference to the
158  global values.
159 
160  Inside the DenseMatrix @a values, the `(i,j)` entry corresponds to the
161  `i`-th vector component at the `j`-th quadrature point.
162  */
163  inline void GetValues(int idx, DenseMatrix &values);
164 
165  /// Return all values associated with mesh element @a idx in a const DenseMatrix.
166  /** The result is stored in the DenseMatrix @a values as a copy of the
167  global values.
168 
169  Inside the DenseMatrix @a values, the `(i,j)` entry corresponds to the
170  `i`-th vector component at the `j`-th quadrature point.
171  */
172  inline void GetValues(int idx, DenseMatrix &values) const;
173 
174  /// Get the IntegrationRule associated with entity (element or face) @a idx.
175  const IntegrationRule &GetIntRule(int idx) const
176  { return GetSpace()->GetIntRule(idx); }
177 
178  /// Write the QuadratureFunction to the stream @a out.
179  void Save(std::ostream &out) const;
180 
181  /// @brief Write the QuadratureFunction to @a out in VTU (ParaView) format.
182  ///
183  /// The data will be uncompressed if @a compression_level is zero, or if the
184  /// format is VTKFormat::ASCII. Otherwise, zlib compression will be used for
185  /// binary data.
186  void SaveVTU(std::ostream &out, VTKFormat format=VTKFormat::ASCII,
187  int compression_level=0) const;
188 
189  /// @brief Save the QuadratureFunction to a VTU (ParaView) file.
190  ///
191  /// The extension ".vtu" will be appended to @a filename.
192  /// @sa SaveVTU(std::ostream &out, VTKFormat format=VTKFormat::ASCII,
193  /// int compression_level=0)
194  void SaveVTU(const std::string &filename, VTKFormat format=VTKFormat::ASCII,
195  int compression_level=0) const;
196 
198  {
199  if (own_qspace) { delete qspace; }
200  }
201 };
202 
203 // Inline methods
204 
206  int idx, Vector &values)
207 {
208  const int s_offset = qspace->offsets[idx];
209  const int sl_size = qspace->offsets[idx+1] - s_offset;
210  values.MakeRef(*this, vdim*s_offset, vdim*sl_size);
211 }
212 
214  int idx, Vector &values) const
215 {
216  const int s_offset = qspace->offsets[idx];
217  const int sl_size = qspace->offsets[idx+1] - s_offset;
218  values.SetSize(vdim*sl_size);
219  values.HostWrite();
220  const double *q = HostRead() + vdim*s_offset;
221  for (int i = 0; i<values.Size(); i++)
222  {
223  values(i) = *(q++);
224  }
225 }
226 
228  int idx, const int ip_num, Vector &values)
229 {
230  const int s_offset = qspace->offsets[idx] * vdim + ip_num * vdim;
231  values.MakeRef(*this, s_offset, vdim);
232 }
233 
235  int idx, const int ip_num, Vector &values) const
236 {
237  const int s_offset = qspace->offsets[idx] * vdim + ip_num * vdim;
238  values.SetSize(vdim);
239  values.HostWrite();
240  const double *q = HostRead() + s_offset;
241  for (int i = 0; i < values.Size(); i++)
242  {
243  values(i) = *(q++);
244  }
245 }
246 
248  int idx, DenseMatrix &values)
249 {
250  const int s_offset = qspace->offsets[idx];
251  const int sl_size = qspace->offsets[idx+1] - s_offset;
252  // Make the values matrix memory an alias of the quadrature function memory
253  Memory<double> &values_mem = values.GetMemory();
254  values_mem.Delete();
255  values_mem.MakeAlias(GetMemory(), vdim*s_offset, vdim*sl_size);
256  values.SetSize(vdim, sl_size);
257 }
258 
260  int idx, DenseMatrix &values) const
261 {
262  const int s_offset = qspace->offsets[idx];
263  const int sl_size = qspace->offsets[idx+1] - s_offset;
264  values.SetSize(vdim, sl_size);
265  values.HostWrite();
266  const double *q = HostRead() + vdim*s_offset;
267  for (int j = 0; j<sl_size; j++)
268  {
269  for (int i = 0; i<vdim; i++)
270  {
271  values(i,j) = *(q++);
272  }
273  }
274 }
275 
276 
278  int vdim_)
279 {
280  if (qspace_ != qspace)
281  {
282  if (own_qspace) { delete qspace; }
283  qspace = qspace_;
284  own_qspace = false;
285  }
286  vdim = (vdim_ < 0) ? vdim : vdim_;
288 }
289 
291  QuadratureSpaceBase *qspace_, double *qf_data, int vdim_)
292 {
293  if (qspace_ != qspace)
294  {
295  if (own_qspace) { delete qspace; }
296  qspace = qspace_;
297  own_qspace = false;
298  }
299  vdim = (vdim_ < 0) ? vdim : vdim_;
300  NewDataAndSize(qf_data, vdim*qspace->GetSize());
301 }
302 
303 
304 } // namespace mfem
305 
306 #endif
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
void NewDataAndSize(double *d, int s)
Set the Vector data and size, deleting the old data, if owned.
Definition: vector.hpp:162
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
void ProjectGridFunction(const GridFunction &gf)
Evaluate a grid function at each quadrature point.
Definition: qfunction.cpp:56
double * HostWrite()
Shortcut for mfem::Write(GetMemory(), TotalSize(), false).
Definition: densemat.hpp:478
Memory< double > & GetMemory()
Definition: densemat.hpp:117
QuadratureFunction(QuadratureSpaceBase *qspace_, double *qf_data, int vdim_=1)
Create a QuadratureFunction based on the given QuadratureSpaceBase, using the external (host) data...
Definition: qfunction.hpp:57
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:512
void Delete()
Delete the owned pointers and reset the Memory object.
QuadratureSpaceBase * qspace
Associated QuadratureSpaceBase object.
Definition: qfunction.hpp:26
virtual const double * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:452
bool own_qspace
Does this own the associated QuadratureSpaceBase?
Definition: qfunction.hpp:27
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
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:460
int vdim
Vector dimension.
Definition: qfunction.hpp:28
Data arrays will be written in ASCII format.
bool OwnsSpace()
Get the QuadratureSpaceBase ownership flag.
Definition: qfunction.hpp:111
Memory< double > & GetMemory()
Return a reference to the Memory object used by the Vector.
Definition: vector.hpp:230
Vector & operator=(const double *v)
Copy Size() entries from v.
Definition: vector.cpp:124
void SaveVTU(std::ostream &out, VTKFormat format=VTKFormat::ASCII, int compression_level=0) const
Write the QuadratureFunction to out in VTU (ParaView) format.
Definition: qfunction.cpp:116
void SetOwnsSpace(bool own)
Set the QuadratureSpaceBase ownership flag.
Definition: qfunction.hpp:114
QuadratureFunction(QuadratureSpaceBase &qspace_, int vdim_=1)
Create a QuadratureFunction based on the given QuadratureSpaceBase.
Definition: qfunction.hpp:38
int GetVDim() const
Get the vector dimension.
Definition: qfunction.hpp:74
const IntegrationRule & GetIntRule(int idx) const
Get the IntegrationRule associated with entity (element or face) idx.
Definition: qfunction.hpp:175
QuadratureFunction()
Default constructor, results in an empty vector.
Definition: qfunction.hpp:32
VTKFormat
Data array format for VTK and VTU files.
Definition: vtk.hpp:96
void SetVDim(int vdim_)
Set the vector dimension, updating the size by calling Vector::SetSize().
Definition: qfunction.hpp:77
QuadratureFunction & operator=(double value)
Set this equal to a constant value.
Definition: qfunction.cpp:19
Array< int > offsets
Entity quadrature point offset array, of size num_entities + 1.
Definition: qspace.hpp:37
Abstract base class for QuadratureSpace and FaceQuadratureSpace.
Definition: qspace.hpp:24
void Save(std::ostream &out) const
Write the QuadratureFunction to the stream out.
Definition: qfunction.cpp:47
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:72
void GetValues(int idx, Vector &values)
Return all values associated with mesh element idx in a Vector.
Definition: qfunction.hpp:205
void MakeAlias(const Memory &base, int offset, int size)
Create a memory object that points inside the memory object base.
QuadratureSpaceBase * GetSpace()
Get the associated QuadratureSpaceBase object.
Definition: qfunction.hpp:81
void SetSpace(QuadratureSpaceBase *qspace_, int vdim_=-1)
Change the QuadratureSpaceBase and optionally the vector dimension.
Definition: qfunction.hpp:277
QuadratureFunction(const QuadratureFunction &orig)
Copy constructor. The QuadratureSpace ownership flag, own_qspace, in the new object is set to false...
Definition: qfunction.hpp:63
Vector data type.
Definition: vector.hpp:60
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition: vector.hpp:576
QuadratureFunction(QuadratureSpaceBase *qspace_, int vdim_=1)
Create a QuadratureFunction based on the given QuadratureSpaceBase.
Definition: qfunction.hpp:47
int GetSize() const
Return the total number of quadrature points.
Definition: qspace.hpp:54
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
const QuadratureSpaceBase * GetSpace() const
Get the associated QuadratureSpaceBase object (const version).
Definition: qfunction.hpp:84
Represents values or vectors of values at quadrature points on a mesh.
Definition: qfunction.hpp:23