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)
83 M_Assert(modes.size() >= Msize,
84 "The Number of requested modes is larger then the available quantity.");
86 Eigen::MatrixXd M = Eigen::MatrixXd::Zero(Msize, Msize);
87 constexpr bool check_vol = std::is_same<volMesh, GeoMesh>::value
88 || std::is_same<surfaceMesh, GeoMesh>::value;
90 if constexpr(check_vol)
97 for (
int i = 0; i < V.size(); i++)
99 M += V(i) * F.transpose().topRows(Msize).col(i) * F.leftCols(Msize).row(i);
104 M = F.transpose().topRows(Msize) * V.asDiagonal() * F.leftCols(Msize);
107 else if constexpr(std::is_same<pointMesh, GeoMesh>::value)
109 M = F.transpose().topRows(Msize) * F.leftCols(Msize);
112 if (Pstream::parRun())
114 reduce(M, sumOp<Eigen::MatrixXd>());
121 PtrList<GeometricField<scalar, fvPatchField, volMesh >>& modes, label Nmodes,
122 bool consider_volumes);
125 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >>& modes,
127 bool consider_volumes =
false);
130 PtrList<GeometricField<vector, fvPatchField, volMesh >>& modes, label Nmodes,
131 bool consider_volumes);
134 PtrList<GeometricField<tensor, fvPatchField, volMesh >>& modes, label Nmodes,
135 bool consider_volumes);
137template<
class Type,
template<
class>
class PatchField,
class GeoMesh>
139 PtrList<GeometricField<Type, PatchField, GeoMesh >>& modes,
140 PtrList<GeometricField<Type, PatchField, GeoMesh >>& modes2, label Nmodes,
141 bool consider_volumes)
147 Msize = modes.size();
148 Msize2 = modes2.size();
156 M_Assert(modes.size() >= Msize,
157 "The Number of requested modes is larger then the available quantity.");
158 M_Assert(modes2.size() >= Msize2,
159 "The Number of requested modes is larger then the available quantity.");
162 Eigen::MatrixXd M = Eigen::MatrixXd::Zero(Msize, Msize2);
164 if (consider_volumes)
171 for (
int i = 0; i < V.size(); i++)
173 M += V(i) * F.transpose().topRows(Msize).col(i) * F2.leftCols(Msize2).row(i);
179 M = F.transpose().topRows(Msize) * V.asDiagonal() * F2.leftCols(Msize2);
184 M = F.transpose().topRows(Msize) * F2.leftCols(Msize2);
187 if (Pstream::parRun())
189 reduce(M, sumOp<Eigen::MatrixXd>());
196 PtrList<GeometricField<scalar, fvPatchField, volMesh >>& modes,
197 PtrList<GeometricField<scalar, fvPatchField, volMesh >>& modes2, label Nmodes,
198 bool consider_volumes);
201 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >>& modes,
202 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >>& modes2,
204 bool consider_volumes =
false);
207 PtrList<GeometricField<vector, fvPatchField, volMesh >>& modes,
208 PtrList<GeometricField<vector, fvPatchField, volMesh >>& modes2, label Nmodes,
209 bool consider_volumes);
212 PtrList<GeometricField<tensor, fvPatchField, volMesh >>& modes,
213 PtrList<GeometricField<tensor, fvPatchField, volMesh >>& modes2, label Nmodes,
214 bool consider_volumes);
217 PtrList<GeometricField<vector, pointPatchField, pointMesh>>& modes,
219 bool consider_volumes);
221template<
class Type,
template<
class>
class PatchField,
class GeoMesh>
223 PtrList<GeometricField<Type, PatchField, GeoMesh >>
225 PtrList<GeometricField<Type, PatchField, GeoMesh >>& modes2,
226 Eigen::VectorXd weights,
228 bool consider_volumes)
234 Msize = modes.size();
235 Msize2 = modes2.size();
243 M_Assert(modes.size() >= Msize,
244 "The Number of requested modes is larger then the available quantity.");
245 M_Assert(modes2.size() >= Msize2,
246 "The Number of requested modes is larger then the available quantity.");
249 Eigen::MatrixXd M = Eigen::MatrixXd::Zero(Msize, Msize2);
251 if (consider_volumes)
258 for (
int i = 0; i < V.size(); i++)
260 M += (V(i) * weights(i)) * F.transpose().topRows(Msize).col(i) * F2.leftCols(
267 M = F.transpose().topRows(Msize) * V.asDiagonal() * ( weights.asDiagonal()
268 * F2.leftCols(Msize2) );
273 M = F.transpose().topRows(Msize) * weights.asDiagonal() * F2.leftCols(Msize2);
276 if (Pstream::parRun())
278 reduce(M, sumOp<Eigen::MatrixXd>());
285 PtrList<GeometricField<scalar, fvPatchField, volMesh >>& modes,
286 PtrList<GeometricField<scalar, fvPatchField, volMesh >>& modes2,
287 Eigen::VectorXd weights,
289 bool consider_volumes);
292 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >>& modes,
293 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >>& modes2,
294 Eigen::VectorXd weights,
296 bool consider_volumes);
299 PtrList<GeometricField<vector, fvPatchField, volMesh >>& modes,
300 PtrList<GeometricField<vector, fvPatchField, volMesh >>& modes2,
301 Eigen::VectorXd weights,
303 bool consider_volumes);
306 PtrList<GeometricField<tensor, fvPatchField, volMesh >>& modes,
307 PtrList<GeometricField<tensor, fvPatchField, volMesh >>& modes2,
308 Eigen::VectorXd weights,
310 bool consider_volumes);
312template<
class Type,
template<
class>
class PatchField,
class GeoMesh>
314 GeometricField<Type, PatchField, GeoMesh>& snapshot)
316 constexpr bool check_vol = std::is_same<volMesh, GeoMesh>::value
317 || std::is_same<surfaceMesh, GeoMesh>::value;
319 if constexpr(check_vol)
322 label dim = std::nearbyint(snapEigen.rows() / (snapshot.mesh().V()).size());
327 else if constexpr(std::is_same<pointMesh, GeoMesh>::value)
330 label dim = std::nearbyint(snapEigen.rows() / (snapshot.mesh()
331 ().points()).size());
339 GeometricField<scalar, fvPatchField, volMesh>& snapshot);
341 GeometricField<vector, fvPatchField, volMesh>& snapshot);
343 GeometricField<vector, pointPatchField, pointMesh>& snapshot);
346 GeometricField<tensor, fvPatchField, volMesh>& snapshot);
348template<
class Type,
template<
class>
class PatchField,
class GeoMesh>
349Eigen::VectorXd
getCoeffs(GeometricField<Type, PatchField, GeoMesh>&
351 PtrList<GeometricField<Type, PatchField, GeoMesh >>& modes, label Nmodes,
352 bool consider_volumes)
358 Msize = modes.size();
365 M_Assert(modes.size() >= Msize,
366 "The Number of requested modes is larger then the available quantity.");
368 Eigen::MatrixXd M_matrix =
getMassMatrix(modes, Nmodes, consider_volumes);
370 Eigen::VectorXd a(Msize);
371 Eigen::VectorXd b(Msize);
372 constexpr bool check_vol = std::is_same<volMesh, GeoMesh>::value
373 || std::is_same<surfaceMesh, GeoMesh>::value;
375 if constexpr(check_vol)
378 b = F.transpose().topRows(Msize) * V.asDiagonal() * snapEigen;
380 else if constexpr(std::is_same<pointMesh, GeoMesh>::value)
382 b = F.transpose().topRows(Msize) * snapEigen;
385 if (Pstream::parRun())
387 reduce(b, sumOp<Eigen::VectorXd>());
390 a = M_matrix.fullPivLu().solve(b);
395 GeometricField<scalar, fvPatchField, volMesh>&
396 snapshot, PtrList<GeometricField<scalar, fvPatchField, volMesh >>& modes,
398 bool consider_volumes);
401 GeometricField<scalar, fvsPatchField, surfaceMesh>&
402 snapshot, PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >>& modes,
404 bool consider_volumes =
false);
407 GeometricField<vector, fvPatchField, volMesh>&
408 snapshot, PtrList<GeometricField<vector, fvPatchField, volMesh >>& modes,
410 bool consider_volumes);
412 GeometricField<vector, pointPatchField, pointMesh>&
413 snapshot, PtrList<GeometricField<vector, pointPatchField, pointMesh>>& modes,
415 bool consider_volumes);
417template<
class Type,
template<
class>
class PatchField,
class GeoMesh >
418Eigen::MatrixXd
getCoeffs(PtrList<GeometricField<Type, PatchField, GeoMesh >>&
420 PtrList<GeometricField<Type, PatchField, GeoMesh >>& modes, label Nmodes,
421 bool consider_volumes)
427 Msize = modes.size();
434 M_Assert(modes.size() >= Msize,
435 "The Number of requested modes is larger then the available quantity.");
436 Eigen::MatrixXd coeff(Msize, snapshots.size());
438 for (
auto i = 0; i < snapshots.size(); i++)
440 coeff.col(i) =
getCoeffs(snapshots[i], modes, Nmodes, consider_volumes);
447 PtrList<GeometricField<scalar, fvPatchField, volMesh >>&
448 snapshot, PtrList<GeometricField<scalar, fvPatchField, volMesh >>& modes,
450 bool consider_volumes);
453 PtrList<GeometricField<vector, fvPatchField, volMesh >>&
454 snapshot, PtrList<GeometricField<vector, fvPatchField, volMesh >>& modes,
456 bool consider_volumes);
459 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >>&
460 snapshot, PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >>& modes,
462 bool consider_volumes);
464 PtrList<GeometricField<vector, pointPatchField, pointMesh>>&
465 snapshot, PtrList<GeometricField<vector, pointPatchField, pointMesh>>& modes,
467 bool consider_volumes);
473 bool consider_volumes =
true;
474 Eigen::MatrixXd coeffsField =
getCoeffs(field, modes, modes.size(), consider_volumes);
483template volVectorField
project_to_POD_basis(volVectorField& field, PtrList<volVectorField>& modes);
484template volScalarField
project_to_POD_basis(volScalarField& field, PtrList<volScalarField>& modes);
489 T fieldCentered = field;
492 bool consider_volumes =
true;
493 Eigen::MatrixXd coeffsField =
getCoeffs(fieldCentered, modes, modes.size(), consider_volumes);
496 T projField = meanField;
501template volVectorField
project_to_POD_basis(volVectorField& field, PtrList<volVectorField>& modes,
const volVectorField& meanField);
502template volScalarField
project_to_POD_basis(volScalarField& field, PtrList<volScalarField>& modes,
const volScalarField& meanField);
506 acquiredSnapshotsTimes,
507 Eigen::MatrixXd parameters)
509 label parsNum = parameters.cols();
510 label parsSamplesNum = parameters.rows();
511 M_Assert(parsSamplesNum == acquiredSnapshotsTimes.size(),
512 "The list of time instants does not have the same number of vectors as the number of parameters samples");
513 Eigen::MatrixXd comb;
514 label totalSnapshotsNum = 0;
516 for (label k = 0; k < acquiredSnapshotsTimes.size(); k++)
518 totalSnapshotsNum += acquiredSnapshotsTimes[k].size();
521 comb.resize(totalSnapshotsNum, parsNum + 1);
524 for (label j = 0; j < acquiredSnapshotsTimes.size(); j++)
526 for (label k = 0; k < acquiredSnapshotsTimes[j].size(); k++)
528 comb(i, parsNum) = (acquiredSnapshotsTimes[j])(k, 0);
529 comb.block(i, 0, 1, parsNum) = parameters.row(j);
Header file of the ITHACAcoeffsMass file.
static Eigen::MatrixXd field2Eigen(GeometricField< Type, 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).
Eigen::VectorXd repeatElements(Eigen::VectorXd input, int n)
Repeat each element n times after themselves.
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...
void setToZero(T &f1)
Set a field of type vol[Scalar|Vector|Tensor]Field to 0.
T project_to_POD_basis(T &field, PtrList< T > &modes)
Compute the projection of a field in a POD basis.
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...
void addFields(T &f1, const T &f2c, double alpha)
Perform the following operation f1 + f2 * alpha with f1 and f2 being of type vol[Scalar|Vector|Tensor...
Eigen::MatrixXd getMassMatrix(PtrList< GeometricField< Type, PatchField, GeoMesh > > &modes, label Nmodes, bool consider_volumes)
Gets the mass matrix using the eigen routine.
void subtractFields(T &f1, const T &f2)
Perform the substraction (f1 - f2) between two fields of type vol[Scalar|Vector|Tensor]Field and alph...