MFEM  v4.5.2
Finite element discretization library
qfunction.cpp
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 #include "qfunction.hpp"
13 #include "quadinterpolator.hpp"
15 
16 namespace mfem
17 {
18 
20 {
21  Vector::operator=(value);
22  return *this;
23 }
24 
26 {
27  MFEM_ASSERT(qspace && v.Size() == this->Size(), "");
29  return *this;
30 }
31 
32 
34 {
35  const char *msg = "invalid input stream";
36  std::string ident;
37 
38  qspace = new QuadratureSpace(mesh, in);
39  own_qspace = true;
40 
41  in >> ident; MFEM_VERIFY(ident == "VDim:", msg);
42  in >> vdim;
43 
44  Load(in, vdim*qspace->GetSize());
45 }
46 
47 void QuadratureFunction::Save(std::ostream &os) const
48 {
49  GetSpace()->Save(os);
50  os << "VDim: " << vdim << '\n'
51  << '\n';
52  Vector::Print(os, vdim);
53  os.flush();
54 }
55 
57 {
58  SetVDim(gf.VectorDim());
59 
60  if (auto *qs_elem = dynamic_cast<QuadratureSpace*>(qspace))
61  {
62  const FiniteElementSpace &gf_fes = *gf.FESpace();
63  const bool use_tensor_products = UsesTensorBasis(gf_fes);
64  const ElementDofOrdering ordering = use_tensor_products ?
67 
68  // Use element restriction to go from L-vector to E-vector
69  const Operator *R = gf_fes.GetElementRestriction(ordering);
70  Vector e_vec(R->Height());
71  R->Mult(gf, e_vec);
72 
73  // Use quadrature interpolator to go from E-vector to Q-vector
74  const QuadratureInterpolator *qi =
75  gf_fes.GetQuadratureInterpolator(*qs_elem);
77  qi->DisableTensorProducts(!use_tensor_products);
78  qi->Values(e_vec, *this);
79  }
80  else if (auto *qs_face = dynamic_cast<FaceQuadratureSpace*>(qspace))
81  {
82  const FiniteElementSpace &gf_fes = *gf.FESpace();
83  const bool use_tensor_products = UsesTensorBasis(gf_fes);
84  const ElementDofOrdering ordering = use_tensor_products ?
87 
88  const FaceType face_type = qs_face->GetFaceType();
89 
90  // Use element restriction to go from L-vector to E-vector
91  const Operator *R = gf_fes.GetFaceRestriction(
92  ordering, face_type, L2FaceValues::SingleValued);
93  Vector e_vec(R->Height());
94  R->Mult(gf, e_vec);
95 
96  // Use quadrature interpolator to go from E-vector to Q-vector
97  const FaceQuadratureInterpolator *qi =
98  gf_fes.GetFaceQuadratureInterpolator(qspace->GetIntRule(0), face_type);
100  qi->DisableTensorProducts(!use_tensor_products);
101  qi->Values(e_vec, *this);
102  }
103  else
104  {
105  // This branch should be unreachable
106  MFEM_ABORT("Unsupported case.");
107  }
108 }
109 
110 std::ostream &operator<<(std::ostream &os, const QuadratureFunction &qf)
111 {
112  qf.Save(os);
113  return os;
114 }
115 
116 void QuadratureFunction::SaveVTU(std::ostream &os, VTKFormat format,
117  int compression_level) const
118 {
119  os << R"(<VTKFile type="UnstructuredGrid" version="0.1")";
120  if (compression_level != 0)
121  {
122  os << R"( compressor="vtkZLibDataCompressor")";
123  }
124  os << " byte_order=\"" << VTKByteOrder() << "\">\n";
125  os << "<UnstructuredGrid>\n";
126 
127  const char *fmt_str = (format == VTKFormat::ASCII) ? "ascii" : "binary";
128  const char *type_str = (format != VTKFormat::BINARY32) ? "Float64" : "Float32";
129  std::vector<char> buf;
130 
131  Mesh &mesh = *qspace->GetMesh();
132 
133  int np = qspace->GetSize();
134  int ne = mesh.GetNE();
135  int sdim = mesh.SpaceDimension();
136 
137  // For quadrature functions, each point is a vertex cell, so number of cells
138  // is equal to number of points
139  os << "<Piece NumberOfPoints=\"" << np
140  << "\" NumberOfCells=\"" << np << "\">\n";
141 
142  // print out the points
143  os << "<Points>\n";
144  os << "<DataArray type=\"" << type_str
145  << "\" NumberOfComponents=\"3\" format=\"" << fmt_str << "\">\n";
146 
147  Vector pt(sdim);
148  for (int i = 0; i < ne; i++)
149  {
151  const IntegrationRule &ir = GetIntRule(i);
152  for (int j = 0; j < ir.Size(); j++)
153  {
154  T.Transform(ir[j], pt);
155  WriteBinaryOrASCII(os, buf, pt[0], " ", format);
156  if (sdim > 1) { WriteBinaryOrASCII(os, buf, pt[1], " ", format); }
157  else { WriteBinaryOrASCII(os, buf, 0.0, " ", format); }
158  if (sdim > 2) { WriteBinaryOrASCII(os, buf, pt[2], "", format); }
159  else { WriteBinaryOrASCII(os, buf, 0.0, "", format); }
160  if (format == VTKFormat::ASCII) { os << '\n'; }
161  }
162  }
163  if (format != VTKFormat::ASCII)
164  {
165  WriteBase64WithSizeAndClear(os, buf, compression_level);
166  }
167  os << "</DataArray>\n";
168  os << "</Points>\n";
169 
170  // Write cells (each cell is just a vertex)
171  os << "<Cells>\n";
172  // Connectivity
173  os << R"(<DataArray type="Int32" Name="connectivity" format=")"
174  << fmt_str << "\">\n";
175 
176  for (int i=0; i<np; ++i) { WriteBinaryOrASCII(os, buf, i, "\n", format); }
177  if (format != VTKFormat::ASCII)
178  {
179  WriteBase64WithSizeAndClear(os, buf, compression_level);
180  }
181  os << "</DataArray>\n";
182  // Offsets
183  os << R"(<DataArray type="Int32" Name="offsets" format=")"
184  << fmt_str << "\">\n";
185  for (int i=0; i<np; ++i) { WriteBinaryOrASCII(os, buf, i, "\n", format); }
186  if (format != VTKFormat::ASCII)
187  {
188  WriteBase64WithSizeAndClear(os, buf, compression_level);
189  }
190  os << "</DataArray>\n";
191  // Types
192  os << R"(<DataArray type="UInt8" Name="types" format=")"
193  << fmt_str << "\">\n";
194  for (int i = 0; i < np; i++)
195  {
196  uint8_t vtk_cell_type = VTKGeometry::POINT;
197  WriteBinaryOrASCII(os, buf, vtk_cell_type, "\n", format);
198  }
199  if (format != VTKFormat::ASCII)
200  {
201  WriteBase64WithSizeAndClear(os, buf, compression_level);
202  }
203  os << "</DataArray>\n";
204  os << "</Cells>\n";
205 
206  os << "<PointData>\n";
207  os << "<DataArray type=\"" << type_str << "\" Name=\"u\" format=\""
208  << fmt_str << "\" NumberOfComponents=\"" << vdim << "\">\n";
209  for (int i = 0; i < ne; i++)
210  {
211  DenseMatrix vals;
212  GetValues(i, vals);
213  for (int j = 0; j < vals.Size(); ++j)
214  {
215  for (int vd = 0; vd < vdim; ++vd)
216  {
217  WriteBinaryOrASCII(os, buf, vals(vd, j), " ", format);
218  }
219  if (format == VTKFormat::ASCII) { os << '\n'; }
220  }
221  }
222  if (format != VTKFormat::ASCII)
223  {
224  WriteBase64WithSizeAndClear(os, buf, compression_level);
225  }
226  os << "</DataArray>\n";
227  os << "</PointData>\n";
228 
229  os << "</Piece>\n";
230  os << "</UnstructuredGrid>\n";
231  os << "</VTKFile>" << std::endl;
232 }
233 
234 void QuadratureFunction::SaveVTU(const std::string &filename, VTKFormat format,
235  int compression_level) const
236 {
237  std::ofstream f(filename + ".vtu");
238  SaveVTU(f, format, compression_level);
239 }
240 
241 }
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
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
QuadratureSpaceBase * qspace
Associated QuadratureSpaceBase object.
Definition: qfunction.hpp:26
bool own_qspace
Does this own the associated QuadratureSpaceBase?
Definition: qfunction.hpp:27
void Print(std::ostream &out=mfem::out, int width=8) const
Prints vector to stream out.
Definition: vector.cpp:725
void DisableTensorProducts(bool disable=true) const
Disable the use of tensor product evaluations, for tensor-product elements, e.g. quads and hexes...
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 void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
int vdim
Vector dimension.
Definition: qfunction.hpp:28
bool UsesTensorBasis(const FiniteElementSpace &fes)
Return true if the mesh contains only one topology and the elements are tensor elements.
Definition: fespace.hpp:957
Data arrays will be written in ASCII format.
void WriteBinaryOrASCII(std::ostream &os, std::vector< char > &buf, const T &val, const char *suffix, VTKFormat format)
Write either ASCII data to the stream or binary data to the buffer depending on the given format...
Definition: vtk.hpp:145
const ElementRestrictionOperator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
Definition: fespace.cpp:1259
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 Load(std::istream **in, int np, int *dim)
Reads a vector from multiple files.
Definition: vector.cpp:57
FaceType
Definition: mesh.hpp:45
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
Native ordering as defined by the FiniteElement.
int Size() const
For backward compatibility define Size to be synonym of Width()
Definition: densemat.hpp:102
A class that performs interpolation from an E-vector to quadrature point values and/or derivatives (Q...
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
const FaceQuadratureInterpolator * GetFaceQuadratureInterpolator(const IntegrationRule &ir, FaceType type) const
Return a FaceQuadratureInterpolator that interpolates E-vectors to quadrature point values and/or der...
Definition: fespace.cpp:1356
Mesh * GetMesh() const
Returns the mesh.
Definition: qspace.hpp:63
void SetVDim(int vdim_)
Set the vector dimension, updating the size by calling Vector::SetSize().
Definition: qfunction.hpp:77
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:683
A class that performs interpolation from a face E-vector to quadrature point values and/or derivative...
void DisableTensorProducts(bool disable=true) const
Disable the use of tensor product evaluations, for tensor-product elements, e.g. quads and hexes...
QuadratureFunction & operator=(double value)
Set this equal to a constant value.
Definition: qfunction.cpp:19
virtual void Save(std::ostream &out) const =0
Write the QuadratureSpace to the stream out.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
void Values(const Vector &e_vec, Vector &q_val) const
Interpolate the values of the E-vector e_vec at quadrature points.
void Save(std::ostream &out) const
Write the QuadratureFunction to the stream out.
Definition: qfunction.cpp:47
void Values(const Vector &e_vec, Vector &q_val) const
Interpolate the values of the E-vector e_vec at quadrature points.
const IntegrationRule & GetIntRule(int idx) const
Return the IntegrationRule associated with entity idx.
Definition: qspace.hpp:72
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
void GetValues(int idx, Vector &values)
Return all values associated with mesh element idx in a Vector.
Definition: qfunction.hpp:205
int SpaceDimension() const
Definition: mesh.hpp:1048
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:936
QuadratureSpaceBase * GetSpace()
Get the associated QuadratureSpaceBase object.
Definition: qfunction.hpp:81
void SetOutputLayout(QVectorLayout layout) const
Set the desired output Q-vector layout. The default value is QVectorLayout::byNODES.
const QuadratureInterpolator * GetQuadratureInterpolator(const IntegrationRule &ir) const
Return a QuadratureInterpolator that interpolates E-vectors to quadrature point values and/or derivat...
Definition: fespace.cpp:1327
int VectorDim() const
Definition: gridfunc.cpp:324
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:348
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
Definition: fespace.hpp:74
const char * VTKByteOrder()
Determine the byte order and return either "BigEndian" or "LittleEndian".
Definition: vtk.cpp:602
static const int POINT
Definition: vtk.hpp:30
Lexicographic ordering for tensor-product FiniteElements.
std::ostream & operator<<(std::ostream &os, SparseMatrix const &mat)
Definition: sparsemat.hpp:711
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
Vector data type.
Definition: vector.hpp:60
void SetOutputLayout(QVectorLayout layout) const
Set the desired output Q-vector layout. The default value is QVectorLayout::byNODES.
Class representing the storage layout of a QuadratureFunction.
Definition: qspace.hpp:92
int GetSize() const
Return the total number of quadrature points.
Definition: qspace.hpp:54
virtual const FaceRestriction * GetFaceRestriction(ElementDofOrdering e_ordering, FaceType, L2FaceValues mul=L2FaceValues::DoubleValued) const
Return an Operator that converts L-vectors to E-vectors on each face.
Definition: fespace.cpp:1292
Abstract operator.
Definition: operator.hpp:24
void WriteBase64WithSizeAndClear(std::ostream &os, std::vector< char > &buf, int compression_level)
Encode in base 64 (and potentially compress) the given data, write it to the output stream (with a he...
Definition: vtk.cpp:654
Represents values or vectors of values at quadrature points on a mesh.
Definition: qfunction.hpp:23
VDIM x NQPT x NE (values) / VDIM x DIM x NQPT x NE (grads)