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 << std::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(Pstream::myProcNo());
490 int last_s = numberOfFiles(casename,
491 "processor" + name(Pstream::myProcNo()));
492
493 if (index >= last_s - 1)
494 {
495 FatalError
496 << "Error: Index " << index << " is out of range. "
497 << "Maximum available index is " << last_s - 2
498 << exit(FatalError);
499 }
500
501 return GeometricField<Type, PatchField, GeoMesh>
502 (
503 IOobject
504 (
505 field.name(),
506 timename + "/" + name(index + 1),
507 field.mesh(),
508 IOobject::MUST_READ
509 ),
510 field.mesh()
511 );
512 }
513}
514
515template<class Type, template<class> class PatchField, class GeoMesh>
516void read_fields(
517 PtrList<GeometricField<Type, PatchField, GeoMesh >> & Lfield, word Name,
518 fileName casename, int first_snap, int n_snap)
519{
520 ITHACAparameters* para(ITHACAparameters::getInstance());
521 fvMesh& mesh = para->mesh;
522 if (!Pstream::parRun())
523 {
524 Info << "######### Reading the Data for " << Name << " #########" << endl;
525 fileName rootpath(".");
526 int last_s;
527 Foam::Time runTime2(Foam::Time::controlDictName, rootpath, casename);
528 if (first_snap >= runTime2.times().size())
529 {
530 Info << "Error the index of the first snapshot must be smaller than the number of snapshots"
531 << endl;
532 exit(0);
533 }
534
535 if (n_snap == 0)
536 {
537 last_s = runTime2.times().size();
538 }
539 else
540 {
541 last_s = min(runTime2.times().size(), n_snap + 2);
542 }
543 for (int i = 2 + first_snap; i < last_s + first_snap; i++)
544 {
545 GeometricField<Type, PatchField, GeoMesh> tmp_field(
546 IOobject
547 (
548 Name,
549 casename + "/" + runTime2.times()[i].name(),
550 mesh,
551 IOobject::MUST_READ
552 ),
553 mesh
554 );
555 Lfield.append(tmp_field.clone());
556 printProgress(double(i + 1) / (last_s + first_snap));
557 }
558 std::cout << std::endl;
559 }
560 else
561 {
562 Info << "######### Reading the Data for " << Name << " #########" <<
563 endl;
564 word timename(mesh.time().rootPath() + "/" +
565 mesh.time().caseName() );
566 timename = timename.substr(0, timename.find_last_of("\\/"));
567 timename = timename + "/" + casename + "/" + "processor" + name(Pstream::myProcNo());
568 int last_s = numberOfFiles(casename,
569 "processor" + name(Pstream::myProcNo()) + "/");
570
571 if (first_snap > last_s)
572 {
573 Info << "Error the index of the first snapshot must be smaller than the number of snapshots"
574 << endl;
575 exit(0);
576 }
577
578 if (n_snap == 0)
579 {
580 }
581 else
582 {
583 last_s = min(last_s, n_snap + 2);
584 }
585
586 for (int i = 1 + first_snap; i < last_s + first_snap; i++)
587 {
588 GeometricField<Type, PatchField, GeoMesh> tmp_field(
589 IOobject
590 (
591 Name,
592 timename + "/" + name(i),
593 mesh,
594 IOobject::MUST_READ
595 ),
596 mesh
597 );
598 Lfield.append(tmp_field.clone());
599 printProgress(double(i + 1) / (last_s + first_snap));
600 }
601
602 Info << endl;
603 }
604}
605
606template<class Type, template<class> class PatchField, class GeoMesh>
607void read_fields(
608 PtrList<GeometricField<Type, PatchField, GeoMesh >> & Lfield,
609 GeometricField<Type, PatchField, GeoMesh>& field,
610 fileName casename, int first_snap, int n_snap)
611{
612 if (!Pstream::parRun())
613 {
614 Info << "######### Reading the Data for " << field.name() << " #########" <<
615 endl;
616 fileName rootpath(".");
617 Foam::Time runTime2(Foam::Time::controlDictName, rootpath, casename);
618 int last_s;
619
620 if (first_snap >= runTime2.times().size())
621 {
622 Info << "Error the index of the first snapshot must be smaller than the number of snapshots"
623 << endl;
624 exit(0);
625 }
626
627 if (n_snap == 0)
628 {
629 last_s = runTime2.times().size();
630 }
631 else
632 {
633 last_s = min(runTime2.times().size(), n_snap + 2);
634 }
635
636 for (int i = 2 + first_snap; i < last_s + first_snap; i++)
637 {
638 GeometricField<Type, PatchField, GeoMesh> tmp_field(
639 IOobject
640 (
641 field.name(),
642 casename + "/" + runTime2.times()[i].name(),
643 field.mesh(),
644 IOobject::MUST_READ
645 ),
646 field.mesh()
647 );
648 Lfield.append(tmp_field.clone());
649 printProgress(double(i + 1) / (last_s + first_snap));
650 }
651
652 std::cout << std::endl;
653 }
654 else
655 {
656 Info << "######### Reading the Data for " << field.name() << " #########" <<
657 endl;
658 word timename(field.mesh().time().rootPath() + "/" +
659 field.mesh().time().caseName() );
660 timename = timename.substr(0, timename.find_last_of("\\/"));
661 timename = timename + "/" + casename + "/" + "processor" + name(Pstream::myProcNo());
662 int last_s = numberOfFiles(casename,
663 "processor" + name(Pstream::myProcNo()) + "/");
664
665 if (first_snap > last_s)
666 {
667 Info << "Error the index of the first snapshot must be smaller than the number of snapshots"
668 << endl;
669 exit(0);
670 }
671
672 if (n_snap == 0)
673 {
674 }
675 else
676 {
677 last_s = min(last_s, n_snap + 2);
678 }
679
680 for (int i = 1 + first_snap; i < last_s + first_snap; i++)
681 {
682 GeometricField<Type, PatchField, GeoMesh> tmp_field(
683 IOobject
684 (
685 field.name(),
686 timename + "/" + name(i),
687 field.mesh(),
688 IOobject::MUST_READ
689 ),
690 field.mesh()
691 );
692 Lfield.append(tmp_field.clone());
693 printProgress(double(i + 1) / (last_s + first_snap));
694 }
695
696 Info << endl;
697 }
698}
699
700template<class Type, template<class> class PatchField, class GeoMesh>
701void readMiddleFields(
702 PtrList<GeometricField<Type, PatchField, GeoMesh >> & Lfield,
703 GeometricField<Type, PatchField, GeoMesh>& field, fileName casename)
704{
705 int par = 1;
706 M_Assert(ITHACAutilities::check_folder(casename + "/" + name(par)) != 0,
707 "No parameter dependent solutions stored into Offline folder");
708
709 while (ITHACAutilities::check_folder(casename + "/" + name(par)))
710 {
711 read_fields(Lfield, field, casename + "/" + name(par));
712 par++;
713 }
714}
715
716template<class Type, template<class> class PatchField, class GeoMesh>
717void readConvergedFields(
718 PtrList<GeometricField<Type, PatchField, GeoMesh >> & Lfield,
719 GeometricField<Type, PatchField, GeoMesh>& field,
720 fileName casename)
721{
722 int par = 1;
723 M_Assert(ITHACAutilities::check_folder(casename + "/" + name(par)) != 0,
724 "No parameter dependent solutions stored into Offline folder");
725 std::cout << "######### Reading the Data for " << field.name() << " #########"
726 << std::endl;
727
728 while (ITHACAutilities::check_folder(casename + "/" + name(par)))
729 {
730 int last = 1;
731
732 while (ITHACAutilities::check_folder(casename + "/" + name(par) + "/" + name(last)))
733 {
734 last++;
735 }
736
737 GeometricField<Type, PatchField, GeoMesh> tmpField(
738 IOobject
739 (
740 field.name(),
741 casename + "/" + name(par) + "/" + name(last - 1),
742 field.mesh(),
743 IOobject::MUST_READ
744 ),
745 field.mesh()
746 );
747 Lfield.append(tmpField.clone());
748 par++;
749 }
750}
751
752template<class Type, template<class> class PatchField, class GeoMesh>
753void read_last_fields(
754 PtrList<GeometricField<Type, PatchField, GeoMesh >> & Lfield,
755 const GeometricField<Type, PatchField, GeoMesh>& field,
756 const fileName casename)
757{
758 if (!Pstream::parRun())
759 {
760 Info << "######### Reading the Data for " << field.name() << " #########" <<
761 endl;
762 fileName rootpath(".");
763 Foam::Time runTime2(Foam::Time::controlDictName, rootpath, casename);
764 int last_s (runTime2.times().size());
765#if defined(OFVER) && (OFVER >= 2212)
766 Lfield.emplace_back
767 (
768 IOobject
769 (
770 field.name(),
771 casename + "/" + runTime2.times()[last_s - 1].name(),
772 field.mesh(),
773 IOobject::MUST_READ
774 ),
775 field.mesh()
776 );
777#else
778 auto tfld =
779 autoPtr<GeometricField<Type, PatchField, GeoMesh >>::New
780 (
781 IOobject
782 (
783 field.name(),
784 casename + "/" + runTime2.times()[last_s - 1].name(),
785 field.mesh(),
786 IOobject::MUST_READ
787 ),
788 field.mesh()
789 );
790 Lfield.append(std::move(tfld));
791#endif
792 std::cout << std::endl;
793 }
794
795 else
796 {
797 Info << "######### Reading the Data for " << field.name() << " #########" <<
798 endl;
799 word timename(field.mesh().time().rootPath() + "/" +
800 field.mesh().time().caseName() );
801 timename = timename.substr(0, timename.find_last_of("\\/"));
802 timename = timename + "/" + casename + "/" + "processor" + name(Pstream::myProcNo());
803 int last_s = numberOfFiles(casename,
804 "processor" + name(Pstream::myProcNo()));
805#if defined(OFVER) && (OFVER >= 2212)
806 Lfield.emplace_back
807 (
808 IOobject
809 (
810 field.name(),
811 timename + "/" + name(last_s - 1),
812 field.mesh(),
813 IOobject::MUST_READ
814 ),
815 field.mesh()
816 );
817#else
818 auto tfld =
819 autoPtr<GeometricField<Type, PatchField, GeoMesh >>::New
820 (
821 IOobject
822 (
823 field.name(),
824 timename + "/" + name(last_s - 1),
825 field.mesh(),
826 IOobject::MUST_READ
827 ),
828 field.mesh()
829 );
830 Lfield.append(std::move(tfld));
831#endif
832 Info << endl;
833 }
834}
835
836template<class Type, template<class> class PatchField, class GeoMesh>
837void readLastFields(
838 PtrList<GeometricField<Type, PatchField, GeoMesh >> & Lfield,
839 const GeometricField<Type, PatchField, GeoMesh>& field,
840 const fileName casename)
841{
842 int par = 1;
843 M_Assert(ITHACAutilities::check_folder(casename + "/" + name(par)) != 0,
844 "No parameter dependent solutions stored into Offline folder");
845
846 while (ITHACAutilities::check_folder(casename + "/" + name(par)))
847 {
848 read_last_fields(Lfield, field, casename + "/" + name(par));
849 par++;
850 }
851}
852
853int numberOfFiles(word folder, word MatrixName, word ext)
854{
855 int number_of_files = 0;
856 std::ifstream in;
857 in.open(folder + "/" + MatrixName + name(0) + ext,
858 std::ios::in | std::ios::binary);
859 Info << folder + "/" + MatrixName + name(0) + ext << endl;
860
861 while (in.good())
862 {
863 in.close();
864 number_of_files++;
865 in.open(folder + "/" + MatrixName + name(number_of_files) + ext,
866 std::ios::in | std::ios::binary);
867 }
868
869 in.close();
870 return number_of_files;
871}
872
873template<class Type, template<class> class PatchField, class GeoMesh>
874void exportFields(
875 PtrList<GeometricField<Type, PatchField, GeoMesh >> & field,
876 word folder, word fieldname)
877{
878 ITHACAutilities::createSymLink(folder);
879 Info << "######### Exporting the Data for " << fieldname << " #########" <<
880 endl;
881 for (int j = 0; j < field.size() ; j++)
882 {
883 exportSolution(field[j], name(j + 1), folder, fieldname);
884 printProgress(double(j + 1) / field.size());
885 }
886
887 std::cout << std::endl;
888}
889
890template void exportFields(
891 PtrList<GeometricField<scalar, fvPatchField, volMesh >>& field,
892 word folder, word fieldname);
893template void exportFields(
894 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh >> & field,
895 word folder, word fieldname);
896template void exportFields(
897 PtrList<GeometricField<vector, fvPatchField, volMesh >>& field,
898 word folder, word fieldname);
899template void exportFields(
900 PtrList<GeometricField<tensor, fvPatchField, volMesh >> & field,
901 word folder, word fieldname);
902
903template<class Type, template<class> class PatchField, class GeoMesh >
904void exportSolution(GeometricField<Type, PatchField, GeoMesh>& s,
905 fileName subfolder, fileName folder,
906 word fieldName)
907{
908 if (!Pstream::parRun())
909 {
910 mkDir(folder + "/" + subfolder);
911 ITHACAutilities::createSymLink(folder);
912 GeometricField<Type, PatchField, GeoMesh> act(fieldName, s);
913 fileName fieldname = folder + "/" + subfolder + "/" + fieldName;
914 OFstream os(fieldname);
915 act.writeHeader(os);
916 os << act << endl;
917 }
918 else
919 {
920 mkDir(folder + "/processor" + name(Pstream::myProcNo()) + "/" + subfolder);
921 ITHACAutilities::createSymLink(folder);
922 GeometricField<Type, PatchField, GeoMesh> act(fieldName, s);
923 fileName fieldname = folder + "/processor" + name(Pstream::myProcNo()) + "/" +
924 subfolder + "/" + fieldName;
925 std::cout << fieldname << std::endl;
926 OFstream os(fieldname);
927 act.writeHeader(os);
928 os << act << endl;
929 }
930}
931
932template void exportSolution(
933 GeometricField<scalar, fvPatchField, volMesh>& s,
934 fileName subfolder, fileName folder,
935 word fieldName);
936template void exportSolution(
937 GeometricField<vector, fvPatchField, volMesh>& s,
938 fileName subfolder, fileName folder,
939 word fieldName);
940template void exportSolution(
941 GeometricField<tensor, fvPatchField, volMesh>& s,
942 fileName subfolder, fileName folder,
943 word fieldName);
944template void exportSolution(
945 GeometricField<scalar, fvsPatchField, surfaceMesh>& s,
946 fileName subfolder, fileName folder,
947 word fieldName);
948
949template<class Type, template<class> class PatchField, class GeoMesh>
950void exportSolution(GeometricField<Type, PatchField, GeoMesh>& s,
951 fileName subfolder, fileName folder)
952{
953 if (!Pstream::parRun())
954 {
955 mkDir(folder + "/" + subfolder);
956 ITHACAutilities::createSymLink(folder);
957 fileName fieldname = folder + "/" + subfolder + "/" + s.name();
958 OFstream os(fieldname);
959 s.writeHeader(os);
960 os << s << endl;
961 }
962 else
963 {
964 mkDir(folder + "/processor" + name(Pstream::myProcNo()) + "/" + subfolder);
965 ITHACAutilities::createSymLink(folder);
966 fileName fieldname = folder + "/processor" + name(Pstream::myProcNo()) + "/" +
967 subfolder + "/" + s.name();
968 OFstream os(fieldname);
969 s.writeHeader(os);
970 os << s << endl;
971 }
972}
973
974template void exportSolution(
975 GeometricField<scalar, fvPatchField, volMesh>& s,
976 fileName subfolder, fileName folder);
977template void exportSolution(
978 GeometricField<vector, fvPatchField, volMesh>& s,
979 fileName subfolder, fileName folder);
980template void exportSolution(
981 GeometricField<tensor, fvPatchField, volMesh>& s,
982 fileName subfolder, fileName folder);
983template void exportSolution(
984 GeometricField<scalar, fvsPatchField, surfaceMesh>& s,
985 fileName subfolder, fileName folder);
986
987template void exportSolution(
988 GeometricField<scalar, pointPatchField, pointMesh>& s,
989 fileName subfolder, fileName folder);
990template void exportSolution(
991 GeometricField<vector, pointPatchField, pointMesh>& s,
992 fileName subfolder, fileName folder);
993template void exportSolution(
994 GeometricField<tensor, pointPatchField, pointMesh>& s,
995 fileName subfolder, fileName folder);
996
997void writePoints(pointField points, fileName folder,
998 fileName subfolder)
999{
1000 if (!Pstream::parRun())
1001 {
1002 mkDir(folder + "/" + subfolder);
1003 ITHACAutilities::createSymLink(folder);
1004 fileName fieldname = folder + "/" + subfolder + "/" + "points";
1005 OFstream os(fieldname);
1006 os << "FoamFile \n { \n version 2.0; \n format ascii; \n class vectorField; \n location ""1 / polyMesh""; \n object points; \n }"
1007 << endl;
1008 os << points << endl;
1009 }
1010 else
1011 {
1012 mkDir(folder + "/processor" + name(Pstream::myProcNo()) + "/" + subfolder);
1013 ITHACAutilities::createSymLink(folder);
1014 fileName fieldname = folder + "/processor" + name(Pstream::myProcNo()) + "/" +
1015 subfolder + "/" + "points";
1016 OFstream os(fieldname);
1017 os << "FoamFile \n { \n version 2.0; \n format ascii; \n class vectorField; \n location ""1 / polyMesh""; \n object points; \n }"
1018 << endl;
1019 os << points << endl;
1020 }
1021}
1022
1023void printProgress(double percentage)
1024{
1025 int val = static_cast<int>(percentage * 100);
1026 int lpad = static_cast<int> (percentage * PBWIDTH);
1027 int rpad = PBWIDTH - lpad;
1028
1029 if (Pstream::master())
1030 {
1031 printf ("\r%3d%% [%.*s%*s]", val, lpad, PBSTR, rpad, "");
1032 fflush (stdout);
1033 }
1034}
1035
1036template<typename T>
1037void save(const List<Eigen::SparseMatrix<T >> & MatrixList, word folder,
1038 word MatrixName)
1039{
1040 mkDir(folder);
1041 for (int i = 0; i < MatrixList.size(); i++)
1042 {
1043 word fileName = folder + "/" + MatrixName + name(i) + ".npz";
1044 cnpy::save(MatrixList[i], fileName);
1045 }
1046}
1047
1048template<typename T>
1049void load(List<Eigen::SparseMatrix<T >> & MatrixList, word folder,
1050 word MatrixName)
1051{
1052 int number_of_files = numberOfFiles(folder, MatrixName, ".npz");
1053 std::cout << "Reading the Matrix " + folder + "/" + MatrixName << std::endl;
1054 M_Assert(number_of_files != 0,
1055 "Check if the file you are trying to read exists" );
1056 MatrixList.resize(0);
1057 Eigen::SparseMatrix<T> A;
1058
1059 for (int i = 0; i < number_of_files; i++)
1060 {
1061 cnpy::load(A, folder + "/" + MatrixName + name(i) + ".npz");
1062 MatrixList.append(A);
1063 }
1064}
1065
1066template void read_fields(PtrList<volScalarField>& Lfield,
1067 word Name,
1068 fileName casename, int first_snap, int n_snap);
1069template void read_fields(PtrList<volVectorField>& Lfield,
1070 word Name,
1071 fileName casename, int first_snap, int n_snap);
1072template void read_fields(PtrList<volTensorField>& Lfield,
1073 word Name,
1074 fileName casename, int first_snap, int n_snap);
1075template void read_fields(PtrList<surfaceScalarField>& Lfield,
1076 word Name,
1077 fileName casename, int first_snap, int n_snap);
1078template void read_fields(PtrList<surfaceVectorField>& Lfield,
1079 word Name,
1080 fileName casename, int first_snap, int n_snap);
1081template void read_fields(PtrList<volScalarField>& Lfield,
1082 volScalarField& field, fileName casename, int first_snap, int n_snap);
1083template void read_fields(PtrList<volVectorField>& Lfield,
1084 volVectorField& field, fileName casename, int first_snap, int n_snap);
1085template void read_fields(PtrList<volTensorField>& Lfield,
1086 volTensorField& field, fileName casename, int first_snap, int n_snap);
1087template void read_fields(PtrList<surfaceScalarField>& Lfield,
1088 surfaceScalarField& field, fileName casename, int first_snap, int n_snap);
1089template void read_fields(PtrList<surfaceVectorField>& Lfield,
1090 surfaceVectorField& field, fileName casename, int first_snap, int n_snap);
1091template void readMiddleFields(PtrList<volScalarField>& Lfield,
1092 volScalarField& field, fileName casename);
1093template void readMiddleFields(PtrList<volVectorField>& Lfield,
1094 volVectorField& field, fileName casename);
1095template void readMiddleFields(PtrList<volTensorField>& Lfield,
1096 volTensorField& field, fileName casename);
1097template void readMiddleFields(PtrList<surfaceScalarField>&
1098 Lfield, surfaceScalarField& field, fileName casename);
1099template void readMiddleFields(PtrList<surfaceVectorField>&
1100 Lfield, surfaceVectorField& field, fileName casename);
1101template void readConvergedFields(PtrList<volScalarField>& Lfield,
1102 volScalarField& field, fileName casename);
1103template void readConvergedFields(PtrList<volVectorField>& Lfield,
1104 volVectorField& field, fileName casename);
1105template void readConvergedFields(PtrList<volTensorField>& Lfield,
1106 volTensorField& field, fileName casename);
1107template void readConvergedFields(PtrList<surfaceScalarField>&
1108 Lfield, surfaceScalarField& field, fileName casename);
1109template void readConvergedFields(PtrList<surfaceVectorField>&
1110 Lfield, surfaceVectorField& field, fileName casename);
1111
1112template void read_last_fields(PtrList<volScalarField>& Lfield,
1113 const volScalarField& field, const fileName casename);
1114template void read_last_fields(PtrList<volVectorField>& Lfield,
1115 const volVectorField& field, const fileName casename);
1116template void read_last_fields(PtrList<volTensorField>& Lfield,
1117 const volTensorField& field, const fileName casename);
1118template void read_last_fields(PtrList<surfaceScalarField>& Lfield,
1119 const surfaceScalarField& field, const fileName casename);
1120template void read_last_fields(PtrList<surfaceVectorField>& Lfield,
1121 const surfaceVectorField& field, const fileName casename);
1122template void readLastFields(PtrList<volScalarField>& Lfield,
1123 const volScalarField& field, const fileName casename);
1124template void readLastFields(PtrList<volVectorField>& Lfield,
1125 const volVectorField& field, const fileName casename);
1126template void readLastFields(PtrList<volTensorField>& Lfield,
1127 const volTensorField& field, const fileName casename);
1128template void readLastFields(PtrList<surfaceScalarField>&
1129 Lfield, const surfaceScalarField& field, const fileName casename);
1130template void readLastFields(PtrList<surfaceVectorField>&
1131 Lfield, const surfaceVectorField& field, const fileName casename);
1132
1133template<typename T>
1134void exportList(T& list, word folder, word filename)
1135{
1136 mkDir(folder);
1137 word fieldname = folder + "/" + filename;
1138 OFstream os(fieldname);
1139
1140 for (int i = 0; i < list.size(); i++)
1141 {
1142 os << list[i] << endl;
1143 }
1144}
1145
1146template void exportList(Field<scalar>& list, word folder,
1147 word filename);
1148template void exportList(Field<vector>& list, word folder,
1149 word filename);
1150template void exportList(Field<tensor>& list, word folder,
1151 word filename);
1152
1153template void save(const List<Eigen::SparseMatrix<double >> & MatrixList,
1154 word folder, word MatrixName);
1155
1156template void load(List<Eigen::SparseMatrix<double >>& MatrixList, word folder,
1157 word MatrixName);
1158
1159
1160template GeometricField<scalar, fvPatchField, volMesh>
1161readFieldByIndex(
1162 const GeometricField<scalar, fvPatchField, volMesh>&,
1163 fileName,
1164 label);
1165
1166template GeometricField<vector, fvPatchField, volMesh>
1167readFieldByIndex(
1168 const GeometricField<vector, fvPatchField, volMesh>&,
1169 fileName,
1170 label);
1171
1172template GeometricField<tensor, fvPatchField, volMesh>
1173readFieldByIndex(
1174 const GeometricField<tensor, fvPatchField, volMesh>&,
1175 fileName,
1176 label);
1177
1178template GeometricField<scalar, fvsPatchField, surfaceMesh>
1179readFieldByIndex(
1180 const GeometricField<scalar, fvsPatchField, surfaceMesh>&,
1181 fileName,
1182 label);
1183
1184template GeometricField<vector, fvsPatchField, surfaceMesh>
1185readFieldByIndex(
1186 const GeometricField<vector, fvsPatchField, surfaceMesh>&,
1187 fileName,
1188 label);
1189
1190}
mesh
Foam::fvMesh & mesh
Definition createMesh.H:47
M_Assert
#define M_Assert(Expr, Msg)
Definition ITHACAassert.H:46
PBWIDTH
#define PBWIDTH
Definition ITHACAregularization.H:71
MAXBUFSIZE
#define MAXBUFSIZE
Definition ITHACAregularization.H:69
PBSTR
#define PBSTR
Definition ITHACAregularization.H:70
ITHACAstream.H
Header file of the ITHACAstream class, it contains the implementation of several methods for input ou...
para
ITHACAparameters * para
Definition NLsolve.H:40
Foam2Eigen::fvMatrix2Eigen
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.
ITHACAparameters
Class for the definition of some general parameters, the parameters must be defined from the file ITH...
Definition ITHACAparameters.H:41
ITHACAparameters::getInstance
static ITHACAparameters * getInstance(fvMesh &mesh, Time &localTime)
Gets an instance of ITHACAparameters, to be used if the instance is not existing.
Definition ITHACAparameters.C:45
T
volScalarField & T
Definition createT.H:46
A
volScalarField & A
Definition createConstants.H:144
Eigen::SliceFromTensor
Matrix< VectorType, Dynamic, Dynamic > SliceFromTensor(Eigen::Tensor< VectorType, 3 > &tensor, label dim, label index1)
Definition EigenFunctions.H:476
ITHACAstream
Namespace for input-output manipulation.
Definition ITHACAstream.C:41
ITHACAstream::writePoints
void writePoints(pointField points, fileName folder, fileName subfolder)
Write points of a mesh to a file.
Definition ITHACAstream.C:997
ITHACAstream::SaveSparseMatrix
void SaveSparseMatrix(Eigen::SparseMatrix< T, Nrows, IND > &m, word folder, word MatrixName)
Export an Eigen sparse matrix into bynary format file.
ITHACAstream::readFieldByIndex
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.
Definition ITHACAstream.C:454
ITHACAstream::numberOfFiles
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.
Definition ITHACAstream.C:853
ITHACAstream::readConvergedFields
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.
Definition ITHACAstream.C:717
ITHACAstream::save
void save(const List< Eigen::SparseMatrix< T > > &MatrixList, word folder, word MatrixName)
Definition ITHACAstream.C:1037
ITHACAstream::exportFields
void exportFields(PtrList< GeometricField< Type, PatchField, GeoMesh > > &field, word folder, word fieldname)
Function to export a scalar of vector field.
Definition ITHACAstream.C:874
ITHACAstream::load
void load(List< Eigen::SparseMatrix< T > > &MatrixList, word folder, word MatrixName)
Definition ITHACAstream.C:1049
ITHACAstream::exportSolution
void exportSolution(GeometricField< Type, PatchField, GeoMesh > &s, fileName subfolder, fileName folder, word fieldName)
Export a field to file in a certain folder and subfolder.
Definition ITHACAstream.C:904
ITHACAstream::read_last_fields
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.
Definition ITHACAstream.C:753
ITHACAstream::readMiddleFields
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...
Definition ITHACAstream.C:701
ITHACAstream::exportMatrix
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 ...
Definition ITHACAstream.C:55
ITHACAstream::exportList
void exportList(T &list, word folder, word filename)
Export a list to file.
Definition ITHACAstream.C:1134
ITHACAstream::SaveDenseMatrix
void SaveDenseMatrix(MatrixType &Matrix, word folder, word MatrixName)
Save a dense matrix to a binary format file.
Definition ITHACAstream.H:675
ITHACAstream::exportFvMatrix
void exportFvMatrix(fvMatrix< Type > &Matrix, word folder, word MatrixName)
Export an Fv Matrix to folder together with its source term.
Definition ITHACAstream.C:44
ITHACAstream::exportVector
void exportVector(Eigen::VectorXd &vector, word Name, word type, word folder)
Definition ITHACAstream.C:255
ITHACAstream::readMatrix
List< Eigen::MatrixXd > readMatrix(word folder, word mat_name)
Read a three dimensional matrix from a txt file in Eigen format.
Definition ITHACAstream.C:372
ITHACAstream::read_fields
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.
Definition ITHACAstream.C:516
ITHACAstream::readLastFields
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.
Definition ITHACAstream.C:837
ITHACAstream::exportTensor
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.
Definition ITHACAstream.C:264
ITHACAstream::printProgress
void printProgress(double percentage)
Print progress bar given the percentage.
Definition ITHACAstream.C:1023
ITHACAutilities::check_folder
bool check_folder(word folder)
Checks if a folder exists.
Definition ITHACAsystem.C:89
ITHACAutilities::createSymLink
void createSymLink(word folder)
Creates symbolic links to 0, system and constant.
Definition ITHACAsystem.C:30
cnpy::save
void save(const Eigen::Matrix< T, -1, dim > &mat, const std::string fname)
cnpy::load
Eigen::Matrix< T, -1, dim > load(Eigen::Matrix< T, -1, dim > &mat, const std::string fname)
i
label i
Definition pEqn.H:46
Generated by   doxygen 1.13.2