30template<
class Type,
template<
class>
class PatchField,
class GeoMesh>
32 PtrList<GeometricField<Type, PatchField, GeoMesh>>& modes,
33 Eigen::MatrixXd& coeff_matrix, label Nmodes)
35 PtrList<GeometricField<Type, PatchField, GeoMesh>> rec_field;
38 for (label k = 0; k < coeff_matrix.cols(); k++)
40 for (label
i = 0;
i < Nmodes;
i++)
43 rec_field.append(modes[
i]*coeff_matrix(
i, k));
47 rec_field[k] += modes[
i] * coeff_matrix(
i, k);
54template PtrList<GeometricField<scalar, fvPatchField, volMesh>>
56 PtrList<GeometricField<scalar, fvPatchField, volMesh>>& modes,
57 Eigen::MatrixXd& coeff_matrix, label Nmodes);
58template PtrList<GeometricField<vector, fvPatchField, volMesh>>
60 PtrList<GeometricField<vector, fvPatchField, volMesh>>& modes,
61 Eigen::MatrixXd& coeff_matrix, label Nmodes);
62template PtrList<GeometricField<tensor, fvPatchField, volMesh>>
64 PtrList<GeometricField<tensor, fvPatchField, volMesh>>& modes,
65 Eigen::MatrixXd& coeff_matrix, label Nmodes);
67template<
class Type,
template<
class>
class PatchField,
class GeoMesh>
69 PtrList<GeometricField<Type, PatchField, GeoMesh>>& modes, label Nmodes,
70 bool consider_volumes)
84 "The Number of requested modes is larger then the available quantity.");
91 M = F.transpose().topRows(Msize) * V.asDiagonal() * F.leftCols(Msize);
95 M = F.transpose().topRows(Msize) * F.leftCols(Msize);
98 if (Pstream::parRun())
100 reduce(M, sumOp<Eigen::MatrixXd>());
107 PtrList<GeometricField<scalar, fvPatchField, volMesh>>& modes, label Nmodes,
108 bool consider_volumes);
111 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh>>& modes,
113 bool consider_volumes =
false);
116 PtrList<GeometricField<vector, fvPatchField, volMesh>>& modes, label Nmodes,
117 bool consider_volumes);
119template<
class Type,
template<
class>
class PatchField,
class GeoMesh>
121 GeometricField<Type, PatchField, GeoMesh>& snapshot)
124 label dim = std::nearbyint(snapEigen.rows() / (snapshot.mesh().V()).size());
126 Eigen::VectorXd vol3 = volumes.replicate(dim, 1);
131 GeometricField<scalar, fvPatchField, volMesh>& snapshot);
133 GeometricField<vector, fvPatchField, volMesh>& snapshot);
135template<
class Type,
template<
class>
class PatchField,
class GeoMesh>
136Eigen::VectorXd
getCoeffs(GeometricField<Type, PatchField, GeoMesh>&
138 PtrList<GeometricField<Type, PatchField, GeoMesh>>& modes, label Nmodes,
139 bool consider_volumes)
145 Msize = modes.size();
153 "The Number of requested modes is larger then the available quantity.");
155 Eigen::MatrixXd M_matrix =
getMassMatrix(modes, Nmodes, consider_volumes);
157 Eigen::VectorXd a(Msize);
158 Eigen::VectorXd b(Msize);
160 if (consider_volumes)
163 b = F.transpose().topRows(Msize) * V.asDiagonal() * snapEigen;
167 b = F.transpose().topRows(Msize) * snapEigen;
170 if (Pstream::parRun())
172 reduce(b, sumOp<Eigen::VectorXd>());
175 a = M_matrix.fullPivLu().solve(b);
180 GeometricField<scalar, fvPatchField, volMesh>&
181 snapshot, PtrList<GeometricField<scalar, fvPatchField, volMesh>>& modes,
183 bool consider_volumes);
186 GeometricField<scalar, fvsPatchField, surfaceMesh>&
187 snapshot, PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh>>& modes,
189 bool consider_volumes =
false);
192 GeometricField<vector, fvPatchField, volMesh>&
193 snapshot, PtrList<GeometricField<vector, fvPatchField, volMesh>>& modes,
195 bool consider_volumes);
197template<
class Type,
template<
class>
class PatchField,
class GeoMesh>
198Eigen::MatrixXd
getCoeffs(PtrList<GeometricField<Type, PatchField, GeoMesh>>&
200 PtrList<GeometricField<Type, PatchField, GeoMesh>>& modes, label Nmodes,
201 bool consider_volumes)
207 Msize = modes.size();
215 "The Number of requested modes is larger then the available quantity.");
216 Eigen::MatrixXd coeff(Msize, snapshots.size());
218 for (
auto i = 0;
i < snapshots.size();
i++)
220 coeff.col(
i) =
getCoeffs(snapshots[
i], modes, Nmodes, consider_volumes);
227 PtrList<GeometricField<scalar, fvPatchField, volMesh>>&
228 snapshot, PtrList<GeometricField<scalar, fvPatchField, volMesh>>& modes,
230 bool consider_volumes);
233 PtrList<GeometricField<vector, fvPatchField, volMesh>>&
234 snapshot, PtrList<GeometricField<vector, fvPatchField, volMesh>>& modes,
236 bool consider_volumes);
239 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh>>&
240 snapshot, PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh>>& modes,
242 bool consider_volumes);
245 acquiredSnapshotsTimes,
246 Eigen::MatrixXd parameters)
248 label parsNum = parameters.cols();
249 label parsSamplesNum = parameters.rows();
250 M_Assert(parsSamplesNum == acquiredSnapshotsTimes.size(),
251 "The list of time instants does not have the same number of vectors as the number of parameters samples");
252 Eigen::MatrixXd comb;
253 label totalSnapshotsNum = 0;
255 for (label k = 0; k < acquiredSnapshotsTimes.size(); k++)
257 totalSnapshotsNum += acquiredSnapshotsTimes[k].size();
260 comb.resize(totalSnapshotsNum, parsNum + 1);
263 for (label j = 0; j < acquiredSnapshotsTimes.size(); j++)
265 for (label k = 0; k < acquiredSnapshotsTimes[j].size(); k++)
267 comb(
i, parsNum) = (acquiredSnapshotsTimes[j])(k, 0);
268 comb.block(
i, 0, 1, parsNum) = parameters.row(j);
#define M_Assert(Expr, Msg)
Header file of the ITHACAcoeffsMass file.
static Eigen::VectorXd field2Eigen(GeometricField< tensor, PatchField, GeoMesh > &field)
Convert a vector OpenFOAM field into an Eigen Vector.
static Eigen::MatrixXd PtrList2Eigen(PtrList< GeometricField< Type, PatchField, GeoMesh > > &fields, label Nfields=-1)
Convert a PtrList of snapshots to Eigen matrix (only internal field)
Namespace to implement some useful assign operation of OF fields.
Eigen::VectorXd getCoeffs(GeometricField< Type, PatchField, GeoMesh > &snapshot, PtrList< GeometricField< Type, PatchField, GeoMesh > > &modes, label Nmodes, bool consider_volumes)
Projects a snapshot on a basis function and gets the coefficients of the projection.
PtrList< GeometricField< Type, PatchField, GeoMesh > > reconstructFromCoeff(PtrList< GeometricField< Type, PatchField, GeoMesh > > &modes, Eigen::MatrixXd &coeff_matrix, label Nmodes)
Exact reconstruction using a certain number of modes for a list of fields and the projection coeffici...
Eigen::VectorXd getMassMatrixFV(GeometricField< Type, PatchField, GeoMesh > &snapshot)
Gets a vector containing the volumes of each cell of the mesh.
Eigen::MatrixXd parTimeCombMat(List< Eigen::VectorXd > acquiredSnapshotsTimes, Eigen::MatrixXd parameters)
A method to compute the time-parameter combined matrix whose any single element corresponds to a uniq...
Eigen::MatrixXd getMassMatrix(PtrList< GeometricField< Type, PatchField, GeoMesh > > &modes, label Nmodes, bool consider_volumes)
Gets the mass matrix using the eigen routine.