Loading...
Searching...
No Matches
ITHACAstream.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ██╗████████╗██╗ ██╗ █████╗ ██████╗ █████╗ ███████╗██╗ ██╗
3 ██║╚══██╔══╝██║ ██║██╔══██╗██╔════╝██╔══██╗ ██╔════╝██║ ██║
4 ██║ ██║ ███████║███████║██║ ███████║█████╗█████╗ ██║ ██║
5 ██║ ██║ ██╔══██║██╔══██║██║ ██╔══██║╚════╝██╔══╝ ╚██╗ ██╔╝
6 ██║ ██║ ██║ ██║██║ ██║╚██████╗██║ ██║ ██║ ╚████╔╝
7 ╚═╝ ╚═╝ ╚═╝ ╚═╝╚═╝ ╚═╝ ╚═════╝╚═╝ ╚═╝ ╚═╝ ╚═══╝
8
9 * In real Time Highly Advanced Computational Applications for Finite Volumes
10 * Copyright (C) 2017 by the ITHACA-FV authors
11-------------------------------------------------------------------------------
12
13License
14 This file is part of ITHACA-FV
15
16 ITHACA-FV is free software: you can redistribute it and/or modify
17 it under the terms of the GNU Lesser General Public License as published by
18 the Free Software Foundation, either version 3 of the License, or
19 (at your option) any later version.
20
21 ITHACA-FV is distributed in the hope that it will be useful,
22 but WITHOUT ANY WARRANTY; without even the implied warranty of
23 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 GNU Lesser General Public License for more details.
25
26 You should have received a copy of the GNU Lesser General Public License
27 along with ITHACA-FV. If not, see <http://www.gnu.org/licenses/>.
28
29\*---------------------------------------------------------------------------*/
30
31
32#include "ITHACAstream.H"
33
34
37
38// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
39
40namespace ITHACAstream
41{
42
43template<typename Type>
44void exportFvMatrix(fvMatrix<Type>& Matrix, word folder,
45 word MatrixName)
46{
47 Eigen::SparseMatrix<double> A;
48 Eigen::VectorXd b;
49 Foam2Eigen::fvMatrix2Eigen(Matrix, A, b);
50 SaveSparseMatrix(A, folder + "/", "A_" + MatrixName);
51 SaveDenseMatrix(b, folder + "/", "B_" + MatrixName);
52}
53
54template <typename T, int dim>
55void exportMatrix(Eigen::Matrix < T, -1, dim > & matrix,
56 word Name, word type,
57 word folder)
58{
59 std::string message = "The extension \"" + type +
60 "\" was not implemented. Check the list of possible extensions.";
61 M_Assert(type == "python" || type == "matlab"
62 || type == "eigen", message.c_str()
63 );
64 mkDir(folder);
65 word est;
66
67 if (type == "python")
68 {
69 est = ".py";
70 OFstream str(folder + "/" + Name + "_mat" + est);
71 str << Name << "=np.array([";
72
73 for (int i = 0; i < matrix.rows(); i++)
74 {
75 for (int j = 0; j < matrix.cols(); j++)
76 {
77 if (j == 0)
78 {
79 str << "[" << setprecision(10) << matrix(i, j);
80 }
81 else
82 {
83 str << "," << setprecision(10) << matrix(i, j);
84 }
85 }
86
87 if (i != (matrix.rows() - 1))
88 {
89 str << "]," << endl;
90 }
91 }
92
93 str << "]])" << endl;
94 }
95
96 if (type == "matlab")
97 {
98 est = ".m";
99 OFstream str(folder + "/" + Name + "_mat" + est);
100 str << Name << "=[";
101
102 for (int i = 0; i < matrix.rows(); i++)
103 {
104 for (int j = 0; j < matrix.cols(); j++)
105 {
106 str << " " << setprecision(10) << matrix(i, j);
107 }
108
109 if (i != (matrix.rows() - 1))
110 {
111 str << ";" << endl;
112 }
113 }
114
115 str << "];" << endl;
116 }
117
118 if (type == "eigen")
119 {
120 const static Eigen::IOFormat CSVFormat(6, false, ", ", "\n");
121 std::ofstream ofs;
122 ofs.precision(20);
123 ofs.open (folder + "/" + Name + "_mat.txt");
124
125 for (int i = 0; i < matrix.rows(); i++)
126 {
127 for (int j = 0; j < matrix.cols(); j++)
128 {
129 if (j == 0)
130 {
131 ofs << matrix(i, j);
132 }
133 else
134 {
135 ofs << " " << matrix(i, j);
136 }
137 }
138
139 if (i != (matrix.rows() - 1))
140 {
141 ofs << endl;
142 }
143 }
144
145 ofs.close();
146 }
147}
148
149template void exportMatrix(Eigen::Matrix < double, -1,
150 -1 > & matrix, word Name, word type,
151 word folder);
152
153template void exportMatrix(Eigen::Matrix < int, -1,
154 -1 > & matrix, word Name, word type,
155 word folder);
156
157template void exportMatrix(Eigen::Matrix < float, -1,
158 -1 > & matrix, word Name, word type,
159 word folder);
160
161template void exportMatrix(Eigen::Matrix < double, -1,
162 1 > & matrix, word Name, word type,
163 word folder);
164
165template void exportMatrix(Eigen::Matrix < int, -1,
166 1 > & matrix, word Name, word type,
167 word folder);
168
169template void exportMatrix(Eigen::Matrix < float, -1,
170 1 > & matrix, word Name, word type,
171 word folder);
172
173void exportMatrix(List <Eigen::MatrixXd>& matrix, word Name,
174 word type, word folder)
175{
176 std::string message = "The extension \"" + type +
177 "\" was not implemented. Check the list of possible extensions.";
178 M_Assert(type == "python" || type == "matlab"
179 || type == "eigen", message.c_str()
180 );
181 mkDir(folder);
182 word est;
183
184 // Python Case
185 if (type == "python")
186 {
187 est = ".py";
188 OFstream str(folder + "/" + Name + "_mat" + est);
189 str << Name << "=np.zeros([" << matrix.size() << "," << matrix[0].rows() <<
190 "," << matrix[0].cols() << "])\n";
191
192 for (int i = 0; i < matrix.size(); i++)
193 {
194 str << Name << "[" << i << ",:,:]=np.array([";
195
196 for (int j = 0; j < matrix[0].rows(); j++)
197 {
198 for (int k = 0; k < matrix[0].cols(); k++)
199 {
200 if ( k == 0)
201 {
202 str << "[" << setprecision(10) << matrix[i](j, k);
203 }
204 else
205 {
206 str << "," << setprecision(10) << matrix[i](j, k);
207 }
208 }
209
210 if (j != (matrix[0].rows() - 1))
211 {
212 str << "]," << endl;
213 }
214 }
215
216 str << "]])\n" << endl;
217 }
218 }
219 // Matlab case
220 else if (type == "matlab")
221 {
222 est = ".m";
223 OFstream str(folder + "/" + Name + "_mat" + est);
224
225 for (int i = 0; i < matrix.size(); i++)
226 {
227 str << Name << "(" << i + 1 << ",:,:)=[";
228
229 for (int j = 0; j < matrix[0].rows(); j++)
230 {
231 for (int k = 0; k < matrix[0].cols(); k++)
232 {
233 str << " " << setprecision(10) << matrix[i](j, k);
234 }
235
236 if (j != (matrix[0].rows() - 1))
237 {
238 str << ";" << endl;
239 }
240 }
241
242 str << "];" << endl;
243 }
244 }
245 else if (type == "eigen")
246 {
247 for (int i = 0; i < matrix.size(); i++)
248 {
249 word Namei = Name + name(i);
250 exportMatrix(matrix[i], Namei, "eigen", folder);
251 }
252 }
253}
254
255void exportVector(Eigen::VectorXd& vector,
256 word Name, word type,
257 word folder)
258{
259 Eigen::MatrixXd matrix = vector;
260 exportMatrix(matrix, Name, type, folder);
261}
262
263template<typename T>
264void exportTensor(Eigen::Tensor<T, 3> tensor, word Name,
265 word type, word folder)
266{
267 std::string message = "The extension \"" + type +
268 "\" was not implemented. Check the list of possible extensions.";
269 M_Assert(type == "python" || type == "matlab"
270 || type == "eigen", message.c_str()
271 );
272 mkDir(folder);
273 word est;
274
275 // Python Case
276 if (type == "python")
277 {
278 est = ".py";
279 OFstream str(folder + "/" + Name + "_mat" + est);
280 str << Name << "=np.zeros([" << tensor.dimension(0) << "," <<
281 Eigen::SliceFromTensor(
282 tensor, 0,
283 0).rows() <<
284 "," << Eigen::SliceFromTensor(tensor, 0,
285 0).cols() << "])\n";
286
287 for (int i = 0; i < tensor.dimension(0); i++)
288 {
289 str << Name << "[" << i << ",:,:]=np.array([";
290
291 for (int j = 0; j < Eigen::SliceFromTensor(tensor, 0, 0).rows(); j++)
292 {
293 for (int k = 0; k < Eigen::SliceFromTensor(tensor, 0, 0).cols(); k++)
294 {
295 if ( k == 0)
296 {
297 str << "[" << setprecision(10) << Eigen::SliceFromTensor(tensor, 0,
298 i)(j, k);
299 }
300 else
301 {
302 str << "," << setprecision(10) << Eigen::SliceFromTensor(tensor, 0,
303 i)(j, k);
304 }
305 }
306
307 if (j != (Eigen::SliceFromTensor(tensor, 0,
308 0).rows() - 1))
309 {
310 str << "]," << endl;
311 }
312 }
313
314 str << "]])\n" << endl;
315 }
316 }
317 // Matlab case
318 else if (type == "matlab")
319 {
320 est = ".m";
321 OFstream str(folder + "/" + Name + "_mat" + est);
322
323 for (int i = 0; i < tensor.dimension(0); i++)
324 {
325 str << Name << "(" << i + 1 << ",:,:)=[";
326
327 for (int j = 0; j < Eigen::SliceFromTensor(tensor, 0,
328 0).rows(); j++)
329 {
330 for (int k = 0; k < Eigen::SliceFromTensor(tensor, 0,
331 0).cols(); k++)
332 {
333 str << " " << setprecision(10) << Eigen::SliceFromTensor(tensor, 0,
334 i)(j, k);
335 }
336
337 if (j != (Eigen::SliceFromTensor(tensor, 0,
338 0).rows() - 1))
339 {
340 str << ";" << endl;
341 }
342 }
343
344 str << "];" << endl;
345 }
346 }
347 else if (type == "eigen")
348 {
349 for (int i = 0; i < tensor.dimension(0); i++)
350 {
351 Eigen::Matrix < T, -1, -1 > matrixAux = Eigen::SliceFromTensor(tensor, 0, i);
352 word Namei = Name + name(i);
353 exportMatrix(matrixAux, Namei, "eigen", folder);
354 }
355 }
356}
357
358
359
360template void exportTensor(Eigen::Tensor<double, 3> tensor,
361 word Name,
362 word type, word folder);
363
364template void exportTensor(Eigen::Tensor<int, 3> tensor,
365 word Name,
366 word type, word folder);
367
368template void exportTensor(Eigen::Tensor<float, 3> tensor,
369 word Name,
370 word type, word folder);
371
372List<Eigen::MatrixXd> readMatrix(word folder, word mat_name)
373{
374 int file_count = 0;
375 DIR* dirp;
376 struct dirent* entry;
377 dirp = opendir(folder.c_str());
378 List <Eigen::MatrixXd> result;
379
380 while ((entry = readdir(dirp)) != NULL)
381 {
382 if (entry->d_type == DT_REG)
383 {
384 file_count++;
385 }
386 }
387
388 closedir (dirp);
389
390 for (int j = 0; j < file_count ; j++)
391 {
392 word matname = folder + "/" + mat_name + name(j) + "_mat.txt";
393 Eigen::MatrixXd temp = readMatrix(matname);
394 result.append(temp);
395 }
396
397 return result;
398}
399
400Eigen::MatrixXd readMatrix(word filename)
401{
402 int cols = 0, rows = 0;
403 double buff[MAXBUFSIZE];
404 // Read numbers from file into buffer.
405 std::ifstream infile;
406 infile.open(filename.c_str());
407 std::string message = "The matrix file \"" + filename +
408 "\" does not exist. Check the existence of the file or the way it is named.";
409 M_Assert(infile.good() != 0, message.c_str()
410 );
411
412 while (! infile.eof())
413 {
414 string line;
415 getline(infile, line);
416 int temp_cols = 0;
417 std::stringstream stream(line);
418
419 while (! stream.eof())
420 {
421 stream >> buff[cols * rows + temp_cols++];
422 }
423
424 if (temp_cols == 0)
425 {
426 continue;
427 }
428
429 if (cols == 0)
430 {
431 cols = temp_cols;
432 }
433
434 rows++;
435 }
436
437 infile.close();
438 rows--;
439 // Populate matrix with numbers.
440 Eigen::MatrixXd result(rows, cols);
441
442 for (int i = 0; i < rows; i++)
443 {
444 for (int j = 0; j < cols; j++)
445 {
446 result(i, j) = buff[ cols * i + j ];
447 }
448 }
449
450 return result;
451}
452
453template<class Type, template<class> class PatchField, class GeoMesh>
454GeometricField<Type, PatchField, GeoMesh> readFieldByIndex(
455 const GeometricField<Type, PatchField, GeoMesh>& field,
456 fileName casename,
457 label index)
458{
459 if (!Pstream::parRun())
460 {
461 fileName rootpath(".");
462 Foam::Time runTime2(Foam::Time::controlDictName, rootpath, casename);
463
464 if (index >= runTime2.times().size() - 2)
465 {
466 FatalError
467 << "Error: Index " << index << " is out of range. "
468 << "Maximum available index is " << runTime2.times().size() - 3
469 << exit(FatalError);
470 }
471
472 return GeometricField<Type, PatchField, GeoMesh>
473 (
474 IOobject
475 (
476 field.name(),
477 casename + "/" + runTime2.times()[index + 2].name(),
478 field.mesh(),
479 IOobject::MUST_READ
480 ),
481 field.mesh()
482 );
483 }
484 else
485 {
486 word timename(field.mesh().time().rootPath() + "/" +
487 field.mesh().time().caseName());
488 timename = timename.substr(0, timename.find_last_of("\\/"));
489 timename = timename + "/" + casename + "/" + "processor" + name(
490 Pstream::myProcNo());
491 int last_s = numberOfFiles(casename,
492 "processor" + name(Pstream::myProcNo()));
493
494 if (index >= last_s - 1)
495 {
496 FatalError
497 << "Error: Index " << index << " is out of range. "
498 << "Maximum available index is " << last_s - 2
499 << exit(FatalError);
500 }
501
502 return GeometricField<Type, PatchField, GeoMesh>
503 (
504 IOobject
505 (
506 field.name(),
507 timename + "/" + name(index + 1),
508 field.mesh(),
509 IOobject::MUST_READ
510 ),
511 field.mesh()
512 );
513 }
514}
515
516template<class Type, template<class> class PatchField, class GeoMesh>
518 PtrList<GeometricField<Type, PatchField, GeoMesh >>& Lfield, word Name,
519 fileName casename, int first_snap, int n_snap)
520{
522 const fvMesh& mesh = para->mesh;
523 const pointMesh& pMesh = pointMesh::New(mesh);
524 constexpr bool check_vol = std::is_same<volMesh, GeoMesh>::value
525 || std::is_same<surfaceMesh, GeoMesh>::value;
526
527 if (!Pstream::parRun())
528 {
529 Info << "######### Reading the Data for " << Name << " #########" << endl;
530 fileName rootpath(".");
531 int last_s;
532 Foam::Time runTime2(Foam::Time::controlDictName, rootpath, casename);
533
534 if (first_snap >= runTime2.times().size())
535 {
536 Info << "Error the index of the first snapshot must be smaller than the number of snapshots"
537 << endl;
538 exit(0);
539 }
540
541 if (n_snap == 0)
542 {
543 last_s = runTime2.times().size();
544 }
545 else
546 {
547 last_s = min(runTime2.times().size(), n_snap + 2);
548 }
549
550 if constexpr(check_vol)
551 {
552 for (int i = 2 + first_snap; i < last_s + first_snap; i++)
553 {
554 GeometricField<Type, PatchField, GeoMesh> tmp_field(
555 IOobject
556 (
557 Name,
558 casename + runTime2.times()[i].name(),
559 mesh,
560 IOobject::MUST_READ
561 ),
562 mesh
563 );
564 Lfield.append(tmp_field.clone());
565 printProgress(double(i + 1) / (last_s + first_snap));
566 }
567 }
568 else if constexpr(std::is_same<pointMesh, GeoMesh>::value)
569 {
570 for (int i = 2 + first_snap; i < last_s + first_snap; i++)
571 {
572 GeometricField<Type, PatchField, GeoMesh> tmp_field(
573 IOobject
574 (
575 Name,
576 casename + runTime2.times()[i].name(),
577 mesh,
578 IOobject::MUST_READ
579 ),
580 pMesh
581 );
582 Lfield.append(tmp_field.clone());
583 printProgress(double(i + 1) / (last_s + first_snap));
584 }
585 }
586
587 Info << endl;
588 }
589 else
590 {
591 Info << "######### Reading the Data for " << Name << " #########" <<
592 endl;
593 word timename(mesh.time().rootPath() + "/" +
594 mesh.time().caseName() );
595 timename = timename.substr(0, timename.find_last_of("\\/"));
596 timename = timename + "/" + casename + "/" + "processor" + name(
597 Pstream::myProcNo());
598 int last_s = numberOfFiles(casename,
599 "processor" + name(Pstream::myProcNo()) + "/");
600
601 if (first_snap > last_s)
602 {
603 Info << "Error the index of the first snapshot must be smaller than the number of snapshots"
604 << endl;
605 exit(0);
606 }
607
608 if (n_snap == 0)
609 {
610 }
611 else
612 {
613 last_s = min(last_s, n_snap + 2);
614 }
615
616 if constexpr(check_vol)
617 {
618 for (int i = 1 + first_snap; i < last_s + first_snap; i++)
619 {
620 GeometricField<Type, PatchField, GeoMesh> tmp_field(
621 IOobject
622 (
623 Name,
624 timename + "/" + name(i),
625 mesh,
626 IOobject::MUST_READ
627 ),
628 mesh
629 );
630 Lfield.append(tmp_field.clone());
631 printProgress(double(i + 1) / (last_s + first_snap));
632 }
633 }
634 else if constexpr(std::is_same<pointMesh, GeoMesh>::value)
635 {
636 for (int i = 1 + first_snap; i < last_s + first_snap; i++)
637 {
638 GeometricField<Type, PatchField, GeoMesh> tmp_field(
639 IOobject
640 (
641 Name,
642 timename + "/" + name(i),
643 mesh,
644 IOobject::MUST_READ
645 ),
646 pMesh
647 );
648 Lfield.append(tmp_field.clone());
649 printProgress(double(i + 1) / (last_s + first_snap));
650 }
651 }
652
653 Info << endl;
654 }
655}
656
657template<class Type, template<class> class PatchField, class GeoMesh>
659 PtrList<GeometricField<Type, PatchField, GeoMesh >>& Lfield,
660 GeometricField<Type, PatchField, GeoMesh>& field,
661 fileName casename, int first_snap, int n_snap)
662{
664 const fvMesh& mesh = para->mesh;
665 const pointMesh& pMesh = pointMesh::New(mesh );
666 constexpr bool check_vol = std::is_same<volMesh, GeoMesh>::value
667 || std::is_same<surfaceMesh, GeoMesh>::value;
668
669 if (!Pstream::parRun())
670 {
671 Info << "######### Reading the Data for " << field.name() << " #########" <<
672 endl;
673 fileName rootpath(".");
674 Foam::Time runTime2(Foam::Time::controlDictName, rootpath, casename);
675 int last_s;
676
677 if (first_snap >= runTime2.times().size())
678 {
679 Info << "Error the index of the first snapshot must be smaller than the number of snapshots"
680 << endl;
681 exit(0);
682 }
683
684 if (n_snap == 0)
685 {
686 last_s = runTime2.times().size();
687 }
688 else
689 {
690 last_s = min(runTime2.times().size(), n_snap + 2);
691 }
692
693 if constexpr(check_vol)
694 {
695 for (int i = 2 + first_snap; i < last_s + first_snap; i++)
696 {
697 GeometricField<Type, PatchField, GeoMesh> tmp_field(
698 IOobject
699 (
700 field.name(),
701 casename + runTime2.times()[i].name(),
702 field.mesh(),
703 IOobject::MUST_READ
704 ),
705 field.mesh()
706 );
707 Lfield.append(tmp_field.clone());
708 printProgress(double(i + 1) / (last_s + first_snap));
709 }
710 }
711 else if constexpr(std::is_same<pointMesh, GeoMesh>::value)
712 {
713 for (int i = 2 + first_snap; i < last_s + first_snap; i++)
714 {
715 GeometricField<Type, PatchField, GeoMesh> tmp_field(
716 IOobject
717 (
718 field.name(),
719 casename + runTime2.times()[i].name(),
720 field.mesh().thisDb(),
721 IOobject::MUST_READ
722 ),
723 pMesh
724 );
725 Lfield.append(tmp_field.clone());
726 printProgress(double(i + 1) / (last_s + first_snap));
727 }
728 }
729
730 Info << endl;
731 }
732 else
733 {
734 Info << "################ Parallel Reading the Data for " << field.name() << " #########" <<
735 endl;
736
737 word timename = casename + "processor" + name(Pstream::myProcNo());
738 Foam::Time runTime2(Foam::Time::controlDictName, ".", timename);
739 int last_s = runTime2.times().size();
740
741 timename = field.mesh().time().rootPath() + "/" + field.mesh().time().caseName();
742 timename = timename.substr(0, timename.find_last_of("\\/"));
743 timename = timename + "/" + casename + "/" + "processor" + name(
744 Pstream::myProcNo());
745
746 if (first_snap > last_s)
747 {
748 Info << "Error the index of the first snapshot must be smaller than the number of snapshots"
749 << endl;
750 exit(0);
751 }
752
753 if (n_snap == 0)
754 {
755 }
756 else
757 {
758 last_s = min(last_s, n_snap + 2);
759 }
760
761 if constexpr(check_vol)
762 {
763 for (int i = 2 + first_snap; i < last_s + first_snap; i++)
764 {
765 GeometricField<Type, PatchField, GeoMesh> tmp_field(
766 IOobject
767 (
768 field.name(),
769 timename + "/" + runTime2.times()[i].name(),
770 field.mesh(),
771 IOobject::MUST_READ
772 ),
773 field.mesh()
774 );
775 Lfield.append(tmp_field.clone());
776 printProgress(double(i + 1) / (last_s + first_snap));
777 }
778 }
779
780 Info << endl;
781 }
782}
783
784template<class Type, template<class> class PatchField, class GeoMesh>
786 PtrList<GeometricField<Type, PatchField, GeoMesh >>& Lfield,
787 GeometricField<Type, PatchField, GeoMesh>& field, fileName casename)
788{
789 int par = 1;
790 M_Assert(ITHACAutilities::check_folder(casename + "/" + name(par)) != 0,
791 "No parameter dependent solutions stored into Offline folder");
792
793 while (ITHACAutilities::check_folder(casename + "/" + name(par)))
794 {
795 read_fields(Lfield, field, casename + "/" + name(par));
796 par++;
797 }
798}
799
800template<class Type, template<class> class PatchField, class GeoMesh>
802 PtrList<GeometricField<Type, PatchField, GeoMesh >>& Lfield,
803 GeometricField<Type, PatchField, GeoMesh>& field,
804 fileName casename)
805{
806 int par = 1;
807 M_Assert(ITHACAutilities::check_folder(casename + "/" + name(par)) != 0,
808 "No parameter dependent solutions stored into Offline folder");
809 Info << "######### Reading the Data for " << field.name() << " #########" << endl;
810
811 while (ITHACAutilities::check_folder(casename + "/" + name(par)))
812 {
813 int last = 1;
814
815 while (ITHACAutilities::check_folder(casename + "/" + name(par) + "/" + name(
816 last)))
817 {
818 last++;
819 }
820
821 GeometricField<Type, PatchField, GeoMesh> tmpField(
822 IOobject
823 (
824 field.name(),
825 casename + "/" + name(par) + "/" + name(last - 1),
826 field.mesh(),
827 IOobject::MUST_READ
828 ),
829 field.mesh()
830 );
831 Lfield.append(tmpField.clone());
832 par++;
833 }
834}
835
836template<class Type, template<class> class PatchField, class GeoMesh>
838 PtrList<GeometricField<Type, PatchField, GeoMesh >>& Lfield,
839 const GeometricField<Type, PatchField, GeoMesh>& field,
840 const fileName casename)
841{
842 if (!Pstream::parRun())
843 {
844 Info << "######### Reading the Data for " << field.name() << " #########" <<
845 endl;
846 fileName rootpath(".");
847 Foam::Time runTime2(Foam::Time::controlDictName, rootpath, casename);
848 int last_s (runTime2.times().size());
849#if defined(OFVER) && (OFVER >= 2212)
850 Lfield.emplace_back
851 (
852 IOobject
853 (
854 field.name(),
855 casename + "/" + runTime2.times()[last_s - 1].name(),
856 field.mesh(),
857 IOobject::MUST_READ
858 ),
859 field.mesh()
860 );
861#else
862 auto tfld =
863 autoPtr<GeometricField<Type, PatchField, GeoMesh >>::New
864 (
865 IOobject
866 (
867 field.name(),
868 casename + "/" + runTime2.times()[last_s - 1].name(),
869 field.mesh(),
870 IOobject::MUST_READ
871 ),
872 field.mesh()
873 );
874 Lfield.append(std::move(tfld));
875#endif
876 Info << endl;
877 }
878 else
879 {
880 Info << "######### Reading the Data for " << field.name() << " #########" <<
881 endl;
882 word timename(field.mesh().time().rootPath() + "/" +
883 field.mesh().time().caseName() );
884 timename = timename.substr(0, timename.find_last_of("\\/"));
885 timename = timename + "/" + casename + "/" + "processor" + name(
886 Pstream::myProcNo());
887 int last_s = numberOfFiles(casename,
888 "processor" + name(Pstream::myProcNo()));
889#if defined(OFVER) && (OFVER >= 2212)
890 Lfield.emplace_back
891 (
892 IOobject
893 (
894 field.name(),
895 timename + "/" + name(last_s - 1),
896 field.mesh(),
897 IOobject::MUST_READ
898 ),
899 field.mesh()
900 );
901#else
902 auto tfld =
903 autoPtr<GeometricField<Type, PatchField, GeoMesh >>::New
904 (
905 IOobject
906 (
907 field.name(),
908 timename + "/" + name(last_s - 1),
909 field.mesh(),
910 IOobject::MUST_READ
911 ),
912 field.mesh()
913 );
914 Lfield.append(std::move(tfld));
915#endif
916 Info << endl;
917 }
918}
919
920template<class Type, template<class> class PatchField, class GeoMesh>
922 PtrList<GeometricField<Type, PatchField, GeoMesh >>& Lfield,
923 const GeometricField<Type, PatchField, GeoMesh>& field,
924 const fileName casename)
925{
926 int par = 1;
927 M_Assert(ITHACAutilities::check_folder(casename + "/" + name(par)) != 0,
928 "No parameter dependent solutions stored into Offline folder");
929
930 while (ITHACAutilities::check_folder(casename + "/" + name(par)))
931 {
932 read_last_fields(Lfield, field, casename + "/" + name(par));
933 par++;
934 }
935}
936
937int numberOfFiles(word folder, word MatrixName, word ext)
938{
939 int number_of_files = 0;
940 std::ifstream in;
941 in.open(folder + "/" + MatrixName + name(0) + ext,
942 std::ios::in | std::ios::binary);
943 Info << folder + "/" + MatrixName + name(0) + ext << endl;
944
945 while (in.good())
946 {
947 in.close();
948 number_of_files++;
949 in.open(folder + "/" + MatrixName + name(number_of_files) + ext,
950 std::ios::in | std::ios::binary);
951 }
952
953 in.close();
954 return number_of_files;
955}
956
957template<class Type, template<class> class PatchField, class GeoMesh>
959 PtrList<GeometricField<Type, PatchField, GeoMesh >>& field,
960 word folder, word fieldname)
961{
963 Info << "######### Exporting the Data for " << fieldname << " #########" <<
964 endl;
965
966 for (int j = 0; j < field.size() ; j++)
967 {
968 exportSolution(field[j], name(j + 1), folder, fieldname);
969 printProgress(double(j + 1) / field.size());
970 }
971
972 Info << endl;
973}
974
975template void exportFields(
976 PtrList<GeometricField<scalar, fvPatchField, volMesh >>& field,
977 word folder, word fieldname);
978template void exportFields(
979 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >>& field,
980 word folder, word fieldname);
981template void exportFields(
982 PtrList<GeometricField<vector, fvPatchField, volMesh >>& field,
983 word folder, word fieldname);
984template void exportFields(
985 PtrList<GeometricField<tensor, fvPatchField, volMesh >>& field,
986 word folder, word fieldname);
987template void exportFields(
988 PtrList<GeometricField<vector, pointPatchField, pointMesh>>& field,
989 word folder, word fieldname);
990
991template<class Type, template<class> class PatchField, class GeoMesh >
992void exportSolution(GeometricField<Type, PatchField, GeoMesh>& s,
993 fileName subfolder, fileName folder,
994 word fieldName)
995{
996 if (!Pstream::parRun())
997 {
998 mkDir(folder + "/" + subfolder);
1000 GeometricField<Type, PatchField, GeoMesh> act(fieldName, s);
1001 fileName fieldname = folder + "/" + subfolder + "/" + fieldName;
1002 OFstream os(fieldname);
1003 act.writeHeader(os);
1004 os << act << endl;
1005 }
1006 else
1007 {
1008 mkDir(folder + "/processor" + name(Pstream::myProcNo()) + "/" + subfolder);
1010 GeometricField<Type, PatchField, GeoMesh> act(fieldName, s);
1011 fileName fieldname = folder + "/processor" + name(Pstream::myProcNo()) + "/" +
1012 subfolder + "/" + fieldName;
1013 OFstream os(fieldname);
1014 act.writeHeader(os);
1015 os << act << endl;
1016 }
1017}
1018
1019template void exportSolution(
1020 GeometricField<scalar, fvPatchField, volMesh>& s,
1021 fileName subfolder, fileName folder,
1022 word fieldName);
1023template void exportSolution(
1024 GeometricField<vector, fvPatchField, volMesh>& s,
1025 fileName subfolder, fileName folder,
1026 word fieldName);
1027template void exportSolution(
1028 GeometricField<tensor, fvPatchField, volMesh>& s,
1029 fileName subfolder, fileName folder,
1030 word fieldName);
1031template void exportSolution(
1032 GeometricField<scalar, fvsPatchField, surfaceMesh>& s,
1033 fileName subfolder, fileName folder,
1034 word fieldName);
1035template void exportSolution(
1036 GeometricField<vector, fvsPatchField, surfaceMesh>& s,
1037 fileName subfolder, fileName folder,
1038 word fieldName);
1039
1040template<class Type, template<class> class PatchField, class GeoMesh>
1041void exportSolution(GeometricField<Type, PatchField, GeoMesh>& s,
1042 fileName subfolder, fileName folder)
1043{
1044 if (!Pstream::parRun())
1045 {
1046 mkDir(folder + "/" + subfolder);
1048 fileName fieldname = folder + "/" + subfolder + "/" + s.name();
1049 OFstream os(fieldname);
1050 s.writeHeader(os);
1051 os << s << endl;
1052 }
1053 else
1054 {
1055 mkDir(folder + "/processor" + name(Pstream::myProcNo()) + "/" + subfolder);
1057 fileName fieldname = folder + "/processor" + name(Pstream::myProcNo()) + "/" +
1058 subfolder + "/" + s.name();
1059 OFstream os(fieldname);
1060 s.writeHeader(os);
1061 os << s << endl;
1062 }
1063}
1064
1065template void exportSolution(
1066 GeometricField<scalar, fvPatchField, volMesh>& s,
1067 fileName subfolder, fileName folder);
1068template void exportSolution(
1069 GeometricField<vector, fvPatchField, volMesh>& s,
1070 fileName subfolder, fileName folder);
1071template void exportSolution(
1072 GeometricField<tensor, fvPatchField, volMesh>& s,
1073 fileName subfolder, fileName folder);
1074template void exportSolution(
1075 GeometricField<scalar, fvsPatchField, surfaceMesh>& s,
1076 fileName subfolder, fileName folder);
1077
1078template void exportSolution(
1079 GeometricField<vector, fvsPatchField, surfaceMesh>& s,
1080 fileName subfolder, fileName folder);
1081
1082template void exportSolution(
1083 GeometricField<scalar, pointPatchField, pointMesh>& s,
1084 fileName subfolder, fileName folder);
1085template void exportSolution(
1086 GeometricField<vector, pointPatchField, pointMesh>& s,
1087 fileName subfolder, fileName folder);
1088template void exportSolution(
1089 GeometricField<tensor, pointPatchField, pointMesh>& s,
1090 fileName subfolder, fileName folder);
1091
1092void writePoints(pointField points, fileName folder,
1093 fileName subfolder)
1094{
1095 if (!Pstream::parRun())
1096 {
1097 mkDir(folder + "/" + subfolder);
1099 fileName fieldname = folder + "/" + subfolder + "/" + "points";
1100 OFstream os(fieldname);
1101 os << "FoamFile \n { \n version 2.0; \n format ascii; \n class vectorField; \n location ""1 / polyMesh""; \n object points; \n }"
1102 << endl;
1103 os << points << endl;
1104 }
1105 else
1106 {
1107 mkDir(folder + "/processor" + name(Pstream::myProcNo()) + "/" + subfolder);
1109 fileName fieldname = folder + "/processor" + name(Pstream::myProcNo()) + "/" +
1110 subfolder + "/" + "points";
1111 OFstream os(fieldname);
1112 os << "FoamFile \n { \n version 2.0; \n format ascii; \n class vectorField; \n location ""1 / polyMesh""; \n object points; \n }"
1113 << endl;
1114 os << points << endl;
1115 }
1116}
1117
1118void printProgress(double percentage)
1119{
1120 int val = static_cast<int>(percentage * 100);
1121 int lpad = static_cast<int> (percentage * PBWIDTH);
1122 int rpad = PBWIDTH - lpad;
1123
1124 if (Pstream::master())
1125 {
1126 printf ("\r%3d%% [%.*s%*s]", val, lpad, PBSTR, rpad, "");
1127 fflush (stdout);
1128 }
1129}
1130
1131template<typename T>
1132void save(const List<Eigen::SparseMatrix<T >>& MatrixList, word folder,
1133 word MatrixName)
1134{
1135 mkDir(folder);
1136
1137 for (int i = 0; i < MatrixList.size(); i++)
1138 {
1139 word fileName = folder + "/" + MatrixName + name(i) + ".npz";
1140 cnpy::save(MatrixList[i], fileName);
1141 }
1142}
1143
1144template<typename T>
1145void load(List<Eigen::SparseMatrix<T >>& MatrixList, word folder,
1146 word MatrixName)
1147{
1148 int number_of_files = numberOfFiles(folder, MatrixName, ".npz");
1149 Info << "Reading the Matrix " + folder + "/" + MatrixName << endl;
1150 M_Assert(number_of_files != 0,
1151 "Check if the file you are trying to read exists" );
1152 MatrixList.resize(0);
1153 Eigen::SparseMatrix<T> A;
1154
1155 for (int i = 0; i < number_of_files; i++)
1156 {
1157 cnpy::load(A, folder + "/" + MatrixName + name(i) + ".npz");
1158 MatrixList.append(A);
1159 }
1160}
1161template void read_fields(PtrList<pointVectorField>& Lfield, word Name,
1162 fileName casename, int first_snap, int n_snap);
1163template void read_fields(PtrList<pointVectorField>& Lfield,
1164 pointVectorField& field,
1165 fileName casename, int first_snap, int n_snap);
1166
1167template void read_fields(PtrList<volScalarField>& Lfield,
1168 word Name,
1169 fileName casename, int first_snap, int n_snap);
1170template void read_fields(PtrList<volVectorField>& Lfield,
1171 word Name,
1172 fileName casename, int first_snap, int n_snap);
1173template void read_fields(PtrList<volTensorField>& Lfield,
1174 word Name,
1175 fileName casename, int first_snap, int n_snap);
1176template void read_fields(PtrList<surfaceScalarField>& Lfield,
1177 word Name,
1178 fileName casename, int first_snap, int n_snap);
1179template void read_fields(PtrList<surfaceVectorField>& Lfield,
1180 word Name,
1181 fileName casename, int first_snap, int n_snap);
1182template void read_fields(PtrList<volScalarField>& Lfield,
1183 volScalarField& field, fileName casename, int first_snap, int n_snap);
1184template void read_fields(PtrList<volVectorField>& Lfield,
1185 volVectorField& field, fileName casename, int first_snap, int n_snap);
1186template void read_fields(PtrList<volTensorField>& Lfield,
1187 volTensorField& field, fileName casename, int first_snap, int n_snap);
1188template void read_fields(PtrList<surfaceScalarField>& Lfield,
1189 surfaceScalarField& field, fileName casename, int first_snap, int n_snap);
1190template void read_fields(PtrList<surfaceVectorField>& Lfield,
1191 surfaceVectorField& field, fileName casename, int first_snap, int n_snap);
1192template void readMiddleFields(PtrList<volScalarField>& Lfield,
1193 volScalarField& field, fileName casename);
1194
1195template void readMiddleFields(PtrList<pointVectorField>& Lfield,
1196 pointVectorField& field, fileName casename);
1197template void readMiddleFields(PtrList<volVectorField>& Lfield,
1198 volVectorField& field, fileName casename);
1199template void readMiddleFields(PtrList<volTensorField>& Lfield,
1200 volTensorField& field, fileName casename);
1201template void readMiddleFields(PtrList<surfaceScalarField>&
1202 Lfield, surfaceScalarField& field, fileName casename);
1203template void readMiddleFields(PtrList<surfaceVectorField>&
1204 Lfield, surfaceVectorField& field, fileName casename);
1205template void readConvergedFields(PtrList<volScalarField>& Lfield,
1206 volScalarField& field, fileName casename);
1207template void readConvergedFields(PtrList<volVectorField>& Lfield,
1208 volVectorField& field, fileName casename);
1209template void readConvergedFields(PtrList<volTensorField>& Lfield,
1210 volTensorField& field, fileName casename);
1211template void readConvergedFields(PtrList<surfaceScalarField>&
1212 Lfield, surfaceScalarField& field, fileName casename);
1213template void readConvergedFields(PtrList<surfaceVectorField>&
1214 Lfield, surfaceVectorField& field, fileName casename);
1215
1216template void read_last_fields(PtrList<volScalarField>& Lfield,
1217 const volScalarField& field, const fileName casename);
1218template void read_last_fields(PtrList<volVectorField>& Lfield,
1219 const volVectorField& field, const fileName casename);
1220template void read_last_fields(PtrList<volTensorField>& Lfield,
1221 const volTensorField& field, const fileName casename);
1222template void read_last_fields(PtrList<surfaceScalarField>& Lfield,
1223 const surfaceScalarField& field, const fileName casename);
1224template void read_last_fields(PtrList<surfaceVectorField>& Lfield,
1225 const surfaceVectorField& field, const fileName casename);
1226template void readLastFields(PtrList<volScalarField>& Lfield,
1227 const volScalarField& field, const fileName casename);
1228template void readLastFields(PtrList<volVectorField>& Lfield,
1229 const volVectorField& field, const fileName casename);
1230template void readLastFields(PtrList<volTensorField>& Lfield,
1231 const volTensorField& field, const fileName casename);
1232template void readLastFields(PtrList<surfaceScalarField>&
1233 Lfield, const surfaceScalarField& field, const fileName casename);
1234template void readLastFields(PtrList<surfaceVectorField>&
1235 Lfield, const surfaceVectorField& field, const fileName casename);
1236// template void readLastFields(PtrList<pointVectorField> &
1237// Lfield, const pointVectorField& field, const fileName casename);
1238
1239template<typename T>
1240void exportList(T& list, word folder, word filename)
1241{
1242 mkDir(folder);
1243 word fieldname = folder + "/" + filename;
1244 OFstream os(fieldname);
1245
1246 for (int i = 0; i < list.size(); i++)
1247 {
1248 os << list[i] << endl;
1249 }
1250}
1251
1252template void exportList(Field<scalar>& list, word folder,
1253 word filename);
1254template void exportList(Field<vector>& list, word folder,
1255 word filename);
1256template void exportList(Field<tensor>& list, word folder,
1257 word filename);
1258
1259template void save(const List<Eigen::SparseMatrix<double >>& MatrixList,
1260 word folder, word MatrixName);
1261
1262template void load(List<Eigen::SparseMatrix<double >>& MatrixList, word folder,
1263 word MatrixName);
1264
1265
1266template GeometricField<scalar, fvPatchField, volMesh>
1268 const GeometricField<scalar, fvPatchField, volMesh>&,
1269 fileName,
1270 label);
1271
1272template GeometricField<vector, fvPatchField, volMesh>
1274 const GeometricField<vector, fvPatchField, volMesh>&,
1275 fileName,
1276 label);
1277
1278template GeometricField<tensor, fvPatchField, volMesh>
1280 const GeometricField<tensor, fvPatchField, volMesh>&,
1281 fileName,
1282 label);
1283
1284template GeometricField<scalar, fvsPatchField, surfaceMesh>
1286 const GeometricField<scalar, fvsPatchField, surfaceMesh>&,
1287 fileName,
1288 label);
1289
1290template GeometricField<vector, fvsPatchField, surfaceMesh>
1292 const GeometricField<vector, fvsPatchField, surfaceMesh>&,
1293 fileName,
1294 label);
1295
1296
1297template<typename T>
1298void read_snapshot(T& snapshot, const Foam::word snap_time,
1299 Foam::word path, Foam::word name)
1300{
1301 // ITHACAparameters* para(ITHACAparameters::getInstance());
1302 const fvMesh& mesh = snapshot.mesh();
1303 word arg_path = path;
1304 word pathProcessor = "";
1305
1306 if (name == "default_name")
1307 {
1308 name = snapshot.name();
1309 }
1310
1311 if (Pstream::parRun())
1312 {
1313 // If specific path, changing processor* from casename directory to the target directory
1314 word path_start = mesh.time().rootPath() + "/" + mesh.time().caseName();
1315 path_start = path_start.substr(0, path_start.find_last_of("\\/")) + "/";
1316 pathProcessor = "/processor" + std::to_string(Pstream::myProcNo());
1317
1318 if (!ITHACAutilities::containsSubstring(path, pathProcessor + "/"))
1319 {
1320 path = path_start + arg_path.substr(0, arg_path.find_last_of("\\/")) + pathProcessor
1321 + "/" + arg_path.substr(arg_path.find_last_of("\\/"), arg_path.size());
1322 }
1323 else
1324 {
1325 path = path_start + arg_path;
1326 }
1327 }
1328 else
1329 {
1330 path = arg_path;
1331 }
1332
1333 path = path + "/" + snap_time;
1334
1335 if (!ITHACAutilities::check_file(path + "/" + name))
1336 {
1337 Info << "Error: data not found at :" << endl;
1338 Info << path << endl;
1339 Info << name << endl;
1340 Info << endl;
1341 abort();
1342 }
1343
1344
1345 T snapshot_dummy(
1346 IOobject
1347 (
1348 name,
1349 path,
1350 mesh,
1351 IOobject::MUST_READ
1352 ),
1353 mesh);
1354
1355 snapshot = snapshot_dummy;
1356}
1357
1358
1359template void read_snapshot(Foam::volScalarField& snapshot,
1360 const Foam::word snap_time, Foam::word path, Foam::word name);
1361template void read_snapshot(Foam::volVectorField& snapshot,
1362 const Foam::word snap_time, Foam::word path, Foam::word name);
1363template void read_snapshot(Foam::volTensorField& snapshot,
1364 const Foam::word snap_time, Foam::word path, Foam::word name);
1365
1366
1367
1368}
Header file of the ITHACAstream class, it contains the implementation of several methods for input ou...
static void fvMatrix2Eigen(fvMatrix< type_foam_matrix > foam_matrix, type_A &A, type_B &b)
Convert a FvMatrix OpenFOAM matrix (Linear System) into a Eigen Matrix A and a source vector b.
Class for the definition of some general parameters, the parameters must be defined from the file ITH...
static ITHACAparameters * getInstance(const fvMesh &mesh, Time &localTime)
Gets an instance of ITHACAparameters, to be used if the instance is not existing.
Namespace for input-output manipulation.
void writePoints(pointField points, fileName folder, fileName subfolder)
Write points of a mesh to a file.
void SaveSparseMatrix(Eigen::SparseMatrix< T, Nrows, IND > &m, word folder, word MatrixName)
Export an Eigen sparse matrix into bynary format file.
GeometricField< Type, PatchField, GeoMesh > readFieldByIndex(const GeometricField< Type, PatchField, GeoMesh > &field, fileName casename, label index)
Function to read a single field by index from a folder.
int numberOfFiles(word folder, word MatrixName, word ext)
Count the number of files with a certain name prefix, used by reading functions of list of matrices.
void readConvergedFields(PtrList< GeometricField< Type, PatchField, GeoMesh > > &Lfield, GeometricField< Type, PatchField, GeoMesh > &field, fileName casename)
Function to read a list of volVectorField from name of the field including only converged snapshots.
void exportFields(PtrList< GeometricField< Type, PatchField, GeoMesh > > &field, word folder, word fieldname)
Function to export a scalar of vector field.
void exportSolution(GeometricField< Type, PatchField, GeoMesh > &s, fileName subfolder, fileName folder, word fieldName)
Export a field to file in a certain folder and subfolder.
void read_last_fields(PtrList< GeometricField< Type, PatchField, GeoMesh > > &Lfield, const GeometricField< Type, PatchField, GeoMesh > &field, const fileName casename)
Function to read a list of fields from the name of the field including only the last snapshot.
void readMiddleFields(PtrList< GeometricField< Type, PatchField, GeoMesh > > &Lfield, GeometricField< Type, PatchField, GeoMesh > &field, fileName casename)
Funtion to read a list of volVectorField from name of the field including all the intermediate snapsh...
void exportMatrix(Eigen::Matrix< T, -1, dim > &matrix, word Name, word type, word folder)
Export the reduced matrices in numpy (type=python), matlab (type=matlab) and txt (type=eigen) format ...
void exportList(T &list, word folder, word filename)
Export a list to file.
void SaveDenseMatrix(MatrixType &Matrix, word folder, word MatrixName)
Save a dense matrix to a binary format file.
void exportFvMatrix(fvMatrix< Type > &Matrix, word folder, word MatrixName)
Export an Fv Matrix to folder together with its source term.
List< Eigen::MatrixXd > readMatrix(word folder, word mat_name)
Read a three dimensional matrix from a txt file in Eigen format.
void read_fields(PtrList< GeometricField< Type, PatchField, GeoMesh > > &Lfield, word Name, fileName casename, int first_snap, int n_snap)
Function to read a list of fields from the name of the field and casename.
void readLastFields(PtrList< GeometricField< Type, PatchField, GeoMesh > > &Lfield, const GeometricField< Type, PatchField, GeoMesh > &field, const fileName casename)
Funtion to read a list of volVectorField from name of the field including only the last snapshots.
void exportTensor(Eigen::Tensor< T, 3 > tensor, word Name, word type, word folder)
Export the reduced tensor in numpy (tipo=python), matlab (tipo=matlab) and txt (tipo=eigen) format.
void printProgress(double percentage)
Print progress bar given the percentage.
bool containsSubstring(std::string contain, std::string contained)
Returns true if contained is a substring of contain, false otherwise.
bool check_folder(word folder)
Checks if a folder exists.
bool check_file(std::string fileName)
Function that returns true if a file exists.
void createSymLink(word folder)
Creates symbolic links to 0, system and constant.