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>
454void read_fields(
455 PtrList<GeometricField<Type, PatchField, GeoMesh>>& Lfield, word Name,
456 fileName casename, int first_snap, int n_snap)
457{
458 ITHACAparameters* para(ITHACAparameters::getInstance());
459 fvMesh& mesh = para->mesh;
460
461 if (!Pstream::parRun())
462 {
463 Info << "######### Reading the Data for " << Name << " #########" << endl;
464 fileName rootpath(".");
465 int last_s;
466 Foam::Time runTime2(Foam::Time::controlDictName, rootpath, casename);
467
468 if (first_snap >= runTime2.times().size())
469 {
470 Info << "Error the index of the first snapshot must be smaller than the number of snapshots"
471 << endl;
472 exit(0);
473 }
474
475 if (n_snap == 0)
476 {
477 last_s = runTime2.times().size();
478 }
479 else
480 {
481 last_s = min(runTime2.times().size(), n_snap + 2);
482 }
483
484 for (int i = 2 + first_snap; i < last_s + first_snap; i++)
485 {
486 GeometricField<Type, PatchField, GeoMesh> tmp_field(
487 IOobject
488 (
489 Name,
490 casename + runTime2.times()[i].name(),
491 mesh,
492 IOobject::MUST_READ
493 ),
494 mesh
495 );
496 Lfield.append(tmp_field.clone());
497 printProgress(double(i + 1) / (last_s + first_snap));
498 }
499
500 std::cout << std::endl;
501 }
502 else
503 {
504 Info << "######### Reading the Data for " << Name << " #########" <<
505 endl;
506 word timename(mesh.time().rootPath() + "/" +
507 mesh.time().caseName() );
508 timename = timename.substr(0, timename.find_last_of("\\/"));
509 timename = timename + "/" + casename + "processor" + name(Pstream::myProcNo());
510 int last_s = numberOfFiles(casename,
511 "processor" + name(Pstream::myProcNo()) + "/");
512
513 if (first_snap > last_s)
514 {
515 Info << "Error the index of the first snapshot must be smaller than the number of snapshots"
516 << endl;
517 exit(0);
518 }
519
520 if (n_snap == 0)
521 {
522 }
523 else
524 {
525 last_s = min(last_s, n_snap + 2);
526 }
527
528 for (int i = 1 + first_snap; i < last_s + first_snap; i++)
529 {
530 GeometricField<Type, PatchField, GeoMesh> tmp_field(
531 IOobject
532 (
533 Name,
534 timename + "/" + name(i),
535 mesh,
536 IOobject::MUST_READ
537 ),
538 mesh
539 );
540 Lfield.append(tmp_field.clone());
541 printProgress(double(i + 1) / (last_s + first_snap));
542 }
543
544 Info << endl;
545 }
546}
547
548template<class Type, template<class> class PatchField, class GeoMesh>
549void read_fields(
550 PtrList<GeometricField<Type, PatchField, GeoMesh>>& Lfield,
551 GeometricField<Type, PatchField, GeoMesh>& field,
552 fileName casename, int first_snap, int n_snap)
553{
554 if (!Pstream::parRun())
555 {
556 Info << "######### Reading the Data for " << field.name() << " #########" <<
557 endl;
558 fileName rootpath(".");
559 Foam::Time runTime2(Foam::Time::controlDictName, rootpath, casename);
560 int last_s;
561
562 if (first_snap >= runTime2.times().size())
563 {
564 Info << "Error the index of the first snapshot must be smaller than the number of snapshots"
565 << endl;
566 exit(0);
567 }
568
569 if (n_snap == 0)
570 {
571 last_s = runTime2.times().size();
572 }
573 else
574 {
575 last_s = min(runTime2.times().size(), n_snap + 2);
576 }
577
578 for (int i = 2 + first_snap; i < last_s + first_snap; i++)
579 {
580 GeometricField<Type, PatchField, GeoMesh> tmp_field(
581 IOobject
582 (
583 field.name(),
584 casename + runTime2.times()[i].name(),
585 field.mesh(),
586 IOobject::MUST_READ
587 ),
588 field.mesh()
589 );
590 Lfield.append(tmp_field.clone());
591 printProgress(double(i + 1) / (last_s + first_snap));
592 }
593
594 std::cout << std::endl;
595 }
596 else
597 {
598 Info << "######### Reading the Data for " << field.name() << " #########" <<
599 endl;
600 word timename(field.mesh().time().rootPath() + "/" +
601 field.mesh().time().caseName() );
602 timename = timename.substr(0, timename.find_last_of("\\/"));
603 timename = timename + "/" + casename + "processor" + name(Pstream::myProcNo());
604 int last_s = numberOfFiles(casename,
605 "processor" + name(Pstream::myProcNo()) + "/");
606
607 if (first_snap > last_s)
608 {
609 Info << "Error the index of the first snapshot must be smaller than the number of snapshots"
610 << endl;
611 exit(0);
612 }
613
614 if (n_snap == 0)
615 {
616 }
617 else
618 {
619 last_s = min(last_s, n_snap + 2);
620 }
621
622 for (int i = 1 + first_snap; i < last_s + first_snap; i++)
623 {
624 GeometricField<Type, PatchField, GeoMesh> tmp_field(
625 IOobject
626 (
627 field.name(),
628 timename + "/" + name(i),
629 field.mesh(),
630 IOobject::MUST_READ
631 ),
632 field.mesh()
633 );
634 Lfield.append(tmp_field.clone());
635 printProgress(double(i + 1) / (last_s + first_snap));
636 }
637
638 Info << endl;
639 }
640}
641
642template<class Type, template<class> class PatchField, class GeoMesh>
643void readMiddleFields(
644 PtrList<GeometricField<Type, PatchField, GeoMesh>>& Lfield,
645 GeometricField<Type, PatchField, GeoMesh>& field, fileName casename)
646{
647 int par = 1;
648 M_Assert(ITHACAutilities::check_folder(casename + name(par)) != 0,
649 "No parameter dependent solutions stored into Offline folder");
650
651 while (ITHACAutilities::check_folder(casename + name(par)))
652 {
653 read_fields(Lfield, field, casename + name(par) + "/");
654 par++;
655 }
656}
657
658template<class Type, template<class> class PatchField, class GeoMesh>
659void readConvergedFields(
660 PtrList<GeometricField<Type, PatchField, GeoMesh>>& Lfield,
661 GeometricField<Type, PatchField, GeoMesh>& field,
662 fileName casename)
663{
664 int par = 1;
665 M_Assert(ITHACAutilities::check_folder(casename + name(par)) != 0,
666 "No parameter dependent solutions stored into Offline folder");
667 std::cout << "######### Reading the Data for " << field.name() << " #########"
668 << std::endl;
669
670 while (ITHACAutilities::check_folder(casename + name(par)))
671 {
672 int last = 1;
673
674 while (ITHACAutilities::check_folder(casename + name(par) + "/" + name(last)))
675 {
676 last++;
677 }
678
679 GeometricField<Type, PatchField, GeoMesh> tmpField(
680 IOobject
681 (
682 field.name(),
683 casename + name(par) + "/" + name(last - 1),
684 field.mesh(),
685 IOobject::MUST_READ
686 ),
687 field.mesh()
688 );
689 Lfield.append(tmpField.clone());
690 par++;
691 }
692}
693
694int numberOfFiles(word folder, word MatrixName, word ext)
695{
696 int number_of_files = 0;
697 std::ifstream in;
698 in.open(folder + "/" + MatrixName + name(0) + ext,
699 std::ios::in | std::ios::binary);
700 Info << folder + "/" + MatrixName + name(0) + ext << endl;
701
702 while (in.good())
703 {
704 in.close();
705 number_of_files++;
706 in.open(folder + "/" + MatrixName + name(number_of_files) + ext,
707 std::ios::in | std::ios::binary);
708 }
709
710 in.close();
711 return number_of_files;
712}
713
714template<class Type, template<class> class PatchField, class GeoMesh>
715void exportFields(
716 PtrList<GeometricField<Type, PatchField, GeoMesh>>& field,
717 word folder, word fieldname)
718{
719 ITHACAutilities::createSymLink(folder);
720 Info << "######### Exporting the Data for " << fieldname << " #########" <<
721 endl;
722
723 for (int j = 0; j < field.size() ; j++)
724 {
725 exportSolution(field[j], name(j + 1), folder, fieldname);
726 printProgress(double(j + 1) / field.size());
727 }
728
729 std::cout << std::endl;
730}
731
732template void exportFields(
733 PtrList<GeometricField<scalar, fvPatchField, volMesh>>& field,
734 word folder, word fieldname);
735template void exportFields(
736 PtrList<GeometricField<scalar, fvsPatchField, surfaceMesh>>& field,
737 word folder, word fieldname);
738template void exportFields(
739 PtrList<GeometricField<vector, fvPatchField, volMesh>>& field,
740 word folder, word fieldname);
741template void exportFields(
742 PtrList<GeometricField<tensor, fvPatchField, volMesh>>& field,
743 word folder, word fieldname);
744
745template<class Type, template<class> class PatchField, class GeoMesh>
746void exportSolution(GeometricField<Type, PatchField, GeoMesh>& s,
747 fileName subfolder, fileName folder,
748 word fieldName)
749{
750 if (!Pstream::parRun())
751 {
752 mkDir(folder + "/" + subfolder);
753 ITHACAutilities::createSymLink(folder);
754 GeometricField<Type, PatchField, GeoMesh> act(fieldName, s);
755 fileName fieldname = folder + "/" + subfolder + "/" + fieldName;
756 OFstream os(fieldname);
757 act.writeHeader(os);
758 os << act << endl;
759 }
760 else
761 {
762 mkDir(folder + "/processor" + name(Pstream::myProcNo()) + "/" + subfolder);
763 ITHACAutilities::createSymLink(folder);
764 GeometricField<Type, PatchField, GeoMesh> act(fieldName, s);
765 fileName fieldname = folder + "/processor" + name(Pstream::myProcNo()) + "/" +
766 subfolder + "/" + fieldName;
767 std::cout << fieldname << std::endl;
768 OFstream os(fieldname);
769 act.writeHeader(os);
770 os << act << endl;
771 }
772}
773
774template void exportSolution(
775 GeometricField<scalar, fvPatchField, volMesh>& s,
776 fileName subfolder, fileName folder,
777 word fieldName);
778template void exportSolution(
779 GeometricField<vector, fvPatchField, volMesh>& s,
780 fileName subfolder, fileName folder,
781 word fieldName);
782template void exportSolution(
783 GeometricField<tensor, fvPatchField, volMesh>& s,
784 fileName subfolder, fileName folder,
785 word fieldName);
786template void exportSolution(
787 GeometricField<scalar, fvsPatchField, surfaceMesh>& s,
788 fileName subfolder, fileName folder,
789 word fieldName);
790
791template<class Type, template<class> class PatchField, class GeoMesh>
792void exportSolution(GeometricField<Type, PatchField, GeoMesh>& s,
793 fileName subfolder, fileName folder)
794{
795 if (!Pstream::parRun())
796 {
797 mkDir(folder + "/" + subfolder);
798 ITHACAutilities::createSymLink(folder);
799 fileName fieldname = folder + "/" + subfolder + "/" + s.name();
800 OFstream os(fieldname);
801 s.writeHeader(os);
802 os << s << endl;
803 }
804 else
805 {
806 mkDir(folder + "/processor" + name(Pstream::myProcNo()) + "/" + subfolder);
807 ITHACAutilities::createSymLink(folder);
808 fileName fieldname = folder + "/processor" + name(Pstream::myProcNo()) + "/" +
809 subfolder + "/" + s.name();
810 OFstream os(fieldname);
811 s.writeHeader(os);
812 os << s << endl;
813 }
814}
815
816template void exportSolution(
817 GeometricField<scalar, fvPatchField, volMesh>& s,
818 fileName subfolder, fileName folder);
819template void exportSolution(
820 GeometricField<vector, fvPatchField, volMesh>& s,
821 fileName subfolder, fileName folder);
822template void exportSolution(
823 GeometricField<tensor, fvPatchField, volMesh>& s,
824 fileName subfolder, fileName folder);
825template void exportSolution(
826 GeometricField<scalar, fvsPatchField, surfaceMesh>& s,
827 fileName subfolder, fileName folder);
828
829template void exportSolution(
830 GeometricField<scalar, pointPatchField, pointMesh>& s,
831 fileName subfolder, fileName folder);
832template void exportSolution(
833 GeometricField<vector, pointPatchField, pointMesh>& s,
834 fileName subfolder, fileName folder);
835template void exportSolution(
836 GeometricField<tensor, pointPatchField, pointMesh>& s,
837 fileName subfolder, fileName folder);
838
839void writePoints(pointField points, fileName folder,
840 fileName subfolder)
841{
842 if (!Pstream::parRun())
843 {
844 mkDir(folder + "/" + subfolder);
845 ITHACAutilities::createSymLink(folder);
846 fileName fieldname = folder + "/" + subfolder + "/" + "points";
847 OFstream os(fieldname);
848 os << "FoamFile \n { \n version 2.0; \n format ascii; \n class vectorField; \n location ""1 / polyMesh""; \n object points; \n }"
849 << endl;
850 os << points << endl;
851 }
852 else
853 {
854 mkDir(folder + "/processor" + name(Pstream::myProcNo()) + "/" + subfolder);
855 ITHACAutilities::createSymLink(folder);
856 fileName fieldname = folder + "/processor" + name(Pstream::myProcNo()) + "/" +
857 subfolder + "/" + "points";
858 OFstream os(fieldname);
859 os << "FoamFile \n { \n version 2.0; \n format ascii; \n class vectorField; \n location ""1 / polyMesh""; \n object points; \n }"
860 << endl;
861 os << points << endl;
862 }
863}
864
865void printProgress(double percentage)
866{
867 int val = static_cast<int>(percentage * 100);
868 int lpad = static_cast<int> (percentage * PBWIDTH);
869 int rpad = PBWIDTH - lpad;
870
871 if (Pstream::master())
872 {
873 printf ("\r%3d%% [%.*s%*s]", val, lpad, PBSTR, rpad, "");
874 fflush (stdout);
875 }
876}
877
878template<typename T>
879void save(const List<Eigen::SparseMatrix<T>>& MatrixList, word folder,
880 word MatrixName)
881{
882 mkDir(folder);
883
884 for (int i = 0; i < MatrixList.size(); i++)
885 {
886 word fileName = folder + "/" + MatrixName + name(i) + ".npz";
887 cnpy::save(MatrixList[i], fileName);
888 }
889}
890
891template<typename T>
892void load(List<Eigen::SparseMatrix<T>>& MatrixList, word folder,
893 word MatrixName)
894{
895 int number_of_files = numberOfFiles(folder, MatrixName, ".npz");
896 std::cout << "Reading the Matrix " + folder + "/" + MatrixName << std::endl;
897 M_Assert(number_of_files != 0,
898 "Check if the file you are trying to read exists" );
899 MatrixList.resize(0);
900 Eigen::SparseMatrix<T> A;
901
902 for (int i = 0; i < number_of_files; i++)
903 {
904 cnpy::load(A, folder + "/" + MatrixName + name(i) + ".npz");
905 MatrixList.append(A);
906 }
907}
908
909template void read_fields(PtrList<volScalarField>& Lfield,
910 word Name,
911 fileName casename, int first_snap, int n_snap);
912template void read_fields(PtrList<volVectorField>& Lfield,
913 word Name,
914 fileName casename, int first_snap, int n_snap);
915template void read_fields(PtrList<volTensorField>& Lfield,
916 word Name,
917 fileName casename, int first_snap, int n_snap);
918template void read_fields(PtrList<surfaceScalarField>& Lfield,
919 word Name,
920 fileName casename, int first_snap, int n_snap);
921template void read_fields(PtrList<surfaceVectorField>& Lfield,
922 word Name,
923 fileName casename, int first_snap, int n_snap);
924template void read_fields(PtrList<volScalarField>& Lfield,
925 volScalarField& field, fileName casename, int first_snap, int n_snap);
926template void read_fields(PtrList<volVectorField>& Lfield,
927 volVectorField& field, fileName casename, int first_snap, int n_snap);
928template void read_fields(PtrList<volTensorField>& Lfield,
929 volTensorField& field, fileName casename, int first_snap, int n_snap);
930template void read_fields(PtrList<surfaceScalarField>& Lfield,
931 surfaceScalarField& field, fileName casename, int first_snap, int n_snap);
932template void read_fields(PtrList<surfaceVectorField>& Lfield,
933 surfaceVectorField& field, fileName casename, int first_snap, int n_snap);
934template void readMiddleFields(PtrList<volScalarField>& Lfield,
935 volScalarField& field, fileName casename);
936template void readMiddleFields(PtrList<volVectorField>& Lfield,
937 volVectorField& field, fileName casename);
938template void readMiddleFields(PtrList<volTensorField>& Lfield,
939 volTensorField& field, fileName casename);
940template void readMiddleFields(PtrList<surfaceScalarField>&
941 Lfield, surfaceScalarField& field, fileName casename);
942template void readMiddleFields(PtrList<surfaceVectorField>&
943 Lfield, surfaceVectorField& field, fileName casename);
944template void readConvergedFields(PtrList<volScalarField>& Lfield,
945 volScalarField& field, fileName casename);
946template void readConvergedFields(PtrList<volVectorField>& Lfield,
947 volVectorField& field, fileName casename);
948template void readConvergedFields(PtrList<volTensorField>& Lfield,
949 volTensorField& field, fileName casename);
950template void readConvergedFields(PtrList<surfaceScalarField>&
951 Lfield, surfaceScalarField& field, fileName casename);
952template void readConvergedFields(PtrList<surfaceVectorField>&
953 Lfield, surfaceVectorField& field, fileName casename);
954
955template<typename T>
956void exportList(T& list, word folder, word filename)
957{
958 mkDir(folder);
959 word fieldname = folder + filename;
960 OFstream os(fieldname);
961
962 for (int i = 0; i < list.size(); i++)
963 {
964 os << list[i] << endl;
965 }
966}
967
968template void exportList(Field<scalar>& list, word folder,
969 word filename);
970template void exportList(Field<vector>& list, word folder,
971 word filename);
972template void exportList(Field<tensor>& list, word folder,
973 word filename);
974
975template void save(const List<Eigen::SparseMatrix<double>>& MatrixList,
976 word folder, word MatrixName);
977
978template void load(List<Eigen::SparseMatrix<double>>& MatrixList, word folder,
979 word MatrixName);
980
981}
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...
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::mesh
fvMesh & mesh
Mesh object defined locally.
Definition ITHACAparameters.H:69
ITHACAparameters::getInstance
static ITHACAparameters * getInstance()
Gets an instance of ITHACAparameters, to be used if the instance is already existing.
Definition ITHACAparameters.C:54
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:477
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:839
ITHACAstream::SaveSparseMatrix
void SaveSparseMatrix(Eigen::SparseMatrix< T, Nrows, IND > &m, word folder, word MatrixName)
Export an Eigen sparse matrix into bynary format file.
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:694
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:659
ITHACAstream::save
void save(const List< Eigen::SparseMatrix< T > > &MatrixList, word folder, word MatrixName)
Definition ITHACAstream.C:879
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:715
ITHACAstream::load
void load(List< Eigen::SparseMatrix< T > > &MatrixList, word folder, word MatrixName)
Definition ITHACAstream.C:892
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:746
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:643
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:956
ITHACAstream::SaveDenseMatrix
void SaveDenseMatrix(MatrixType &Matrix, word folder, word MatrixName)
Save a dense matrix to a binary format file.
Definition ITHACAstream.H:617
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:454
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:865
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, std::string order="rowMajor")
i
label i
Definition pEqn.H:46
Generated by   doxygen 1.10.0