Loading...
Searching...
No Matches
Foam2Eigen.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#include "Foam2Eigen.H"
32
35
36// * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * * //
37
38// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
39template <template <class> class PatchField, class GeoMesh>
40Eigen::VectorXd Foam2Eigen::field2Eigen(
41 GeometricField<tensor, PatchField, GeoMesh>& field)
42{
43 Eigen::VectorXd out;
44 out.resize(label(field.size() * 9));
45
46 for (label l = 0; l < field.size(); l++)
47 {
48 for (label j = 0; j < 9; j++)
49 {
50 out(j * field.size() + l) = field[l][j];
51 }
52 }
53
54 return out;
55}
56
57template Eigen::VectorXd Foam2Eigen::field2Eigen(
58 volTensorField& field);
59
60template <template <class> class PatchField, class GeoMesh>
61Eigen::VectorXd Foam2Eigen::field2Eigen(
62 GeometricField<vector, PatchField, GeoMesh>& field)
63{
64 Eigen::VectorXd out;
65 out.resize(label(field.size() * 3));
66
67 for (label l = 0; l < field.size(); l++)
68 {
69 for (label j = 0; j < 3; j++)
70 {
71 out(j * field.size() + l) = field[l][j];
72 }
73 }
74
75 return out;
76}
77
78template Eigen::VectorXd Foam2Eigen::field2Eigen(
79 volVectorField& field);
80
81template <template <class> class PatchField, class GeoMesh>
82Eigen::VectorXd Foam2Eigen::field2Eigen(
83 GeometricField<scalar, PatchField, GeoMesh>& field)
84{
85 Eigen::VectorXd out;
86 out.resize(label(field.size()));
87
88 for (label l = 0; l < field.size(); l++)
89 {
90 out(l) = field[l];
91 }
92
93 return out;
94}
95
96template <template <class> class PatchField, class GeoMesh>
97Eigen::Map<Eigen::MatrixXd> Foam2Eigen::field2EigenMap(
98 GeometricField<scalar, PatchField, GeoMesh>& field)
99{
100 Eigen::Map<Eigen::MatrixXd> output(field.ref().data(), field.size(), 1);
101 return std::move(output);
102}
103
104template <template <class> class PatchField, class GeoMesh>
105Eigen::Map<Eigen::MatrixXd> Foam2Eigen::field2EigenMapBC(
106 GeometricField<scalar, PatchField, GeoMesh>& field, int BC_index)
107{
108 Eigen::Map<Eigen::MatrixXd> output(field.boundaryFieldRef()[BC_index].data(), field.boundaryField()[BC_index].size(), 1);
109 return std::move(output);
110}
111
112template Eigen::Map<Eigen::MatrixXd> Foam2Eigen::field2EigenMap(
113 volScalarField& field);
114
115template Eigen::Map<Eigen::MatrixXd> Foam2Eigen::field2EigenMapBC(
116 volScalarField& field, int BC_index);
117
118template Eigen::VectorXd Foam2Eigen::field2Eigen(
119 volScalarField& field);
120template Eigen::VectorXd Foam2Eigen::field2Eigen(
121 surfaceScalarField& field);
122
123template <>
124Eigen::VectorXd Foam2Eigen::field2Eigen(const Field<scalar>& field)
125{
126 Eigen::VectorXd out;
127 out.resize(label(field.size()));
128
129 for (label l = 0; l < field.size(); l++)
130 {
131 out(l) = field[l];
132 }
133
134 return out;
135}
136
137template <>
138Eigen::VectorXd Foam2Eigen::field2Eigen(const Field<vector>& field)
139{
140 Eigen::VectorXd out;
141 out.resize(label(field.size() * 3));
142
143 for (label l = 0; l < field.size(); l++)
144 {
145 for (label j = 0; j < 3; j++)
146 {
147 out(j * field.size() + l) = field[l][j];
148 }
149 }
150
151 return out;
152}
153
154template <>
155Eigen::VectorXd Foam2Eigen::field2Eigen(const Field<tensor>& field)
156{
157 Eigen::VectorXd out;
158 out.resize(label(field.size() * 9));
159
160 for (label l = 0; l < field.size(); l++)
161 {
162 for (label j = 0; j < 9; j++)
163 {
164 out(j * field.size() + l) = field[l][j];
165 }
166 }
167
168 return out;
169}
170
171template <>
172Eigen::VectorXd Foam2Eigen::field2Eigen(const
173 DimensionedField<scalar, Foam::volMesh>& field)
174{
175 Eigen::VectorXd out;
176 out.resize(label(field.size()));
177
178 for (label l = 0; l < field.size(); l++)
179 {
180 out(l) = field[l];
181 }
182
183 return out;
184}
185
186template <template <class> class PatchField, class GeoMesh>
187List<Eigen::VectorXd> Foam2Eigen::field2EigenBC(
188 GeometricField<tensor, PatchField, GeoMesh>& field)
189{
190 List<Eigen::VectorXd> Out;
191 label size = field.boundaryField().size();
192 Out.resize(size);
193
194 for (label i = 0; i < size; i++)
195 {
196 label sizei = field.boundaryField()[i].size();
197 Out[i].resize(sizei * 9);
198
199 for (label k = 0; k < sizei; k++)
200 {
201 for (label j = 0; j < 9; j++)
202 {
203 Out[i](k + j * sizei) = field.boundaryField()[i][k][j];
204 }
205 }
206 }
207
208 return Out;
209}
210
211template List<Eigen::VectorXd> Foam2Eigen::field2EigenBC(
212 volTensorField& field);
213
214template <template <class> class PatchField, class GeoMesh>
215List<Eigen::VectorXd> Foam2Eigen::field2EigenBC(
216 GeometricField<vector, PatchField, GeoMesh>& field)
217{
218 List<Eigen::VectorXd> Out;
219 label size = field.boundaryField().size();
220 Out.resize(size);
221
222 for (label i = 0; i < size; i++)
223 {
224 label sizei = field.boundaryField()[i].size();
225 Out[i].resize(sizei * 3);
226
227 for (label k = 0; k < sizei; k++)
228 {
229 for (label j = 0; j < 3; j++)
230 {
231 Out[i](k + j * sizei) = field.boundaryField()[i][k][j];
232 }
233 }
234 }
235
236 return Out;
237}
238
239template List<Eigen::VectorXd> Foam2Eigen::field2EigenBC(
240 volVectorField& field);
241
242
243template <template <class> class PatchField, class GeoMesh>
244List<Eigen::VectorXd> Foam2Eigen::field2EigenBC(
245 GeometricField<scalar, PatchField, GeoMesh>& field)
246{
247 List<Eigen::VectorXd> Out;
248 label size = field.boundaryField().size();
249 Out.resize(size);
250
251 for (label i = 0; i < size; i++)
252 {
253 label sizei = field.boundaryField()[i].size();
254 Out[i].resize(sizei);
255
256 for (label k = 0; k < sizei; k++)
257 {
258 Out[i](k) = field.boundaryField()[i][k];
259 }
260 }
261
262 return Out;
263}
264template List<Eigen::VectorXd> Foam2Eigen::field2EigenBC(
265 volScalarField& field);
266
267template <template <class> class PatchField, class GeoMesh>
268List<Eigen::MatrixXd> Foam2Eigen::PtrList2EigenBC(
269 PtrList<GeometricField<scalar, PatchField, GeoMesh>>&
270 fields,
271 label Nfields)
272{
273 label Nf;
274 M_Assert(Nfields <= fields.size(),
275 "The Number of requested fields cannot be bigger than the number of requested entries.");
276
277 if (Nfields == -1)
278 {
279 Nf = fields.size();
280 }
281 else
282 {
283 Nf = Nfields;
284 }
285
286 List<Eigen::MatrixXd> Out;
287 label NBound = fields[0].boundaryField().size();
288 Out.resize(NBound);
289
290 for (label i = 0; i < NBound; i++)
291 {
292 label sizei = fields[0].boundaryField()[i].size();
293 Out[i].resize(sizei, Nf);
294 }
295
296 for (label k = 0; k < Nf; k++)
297 {
298 List<Eigen::VectorXd> temp;
299 temp = field2EigenBC(fields[k]);
300
301 for (label i = 0; i < NBound; i++)
302 {
303 Out[i].col(k) = temp[i];
304 }
305 }
306
307 return Out;
308}
309
310template List<Eigen::MatrixXd> Foam2Eigen::PtrList2EigenBC(
311 PtrList<volScalarField>& fields, label Nfields);
312template List<Eigen::MatrixXd> Foam2Eigen::PtrList2EigenBC(
313 PtrList<surfaceScalarField>& fields, label Nfields);
314
315
316template <template <class> class PatchField, class GeoMesh>
317List<Eigen::MatrixXd> Foam2Eigen::PtrList2EigenBC(
318 PtrList<GeometricField<vector, PatchField, GeoMesh>>&
319 fields,
320 label Nfields)
321{
322 label Nf;
323 M_Assert(Nfields <= fields.size(),
324 "The Number of requested fields cannot be bigger than the number of requested entries.");
325
326 if (Nfields == -1)
327 {
328 Nf = fields.size();
329 }
330 else
331 {
332 Nf = Nfields;
333 }
334
335 List<Eigen::MatrixXd> Out;
336 label NBound = fields[0].boundaryField().size();
337 Out.resize(NBound);
338
339 for (label i = 0; i < NBound; i++)
340 {
341 label sizei = fields[0].boundaryField()[i].size();
342 Out[i].resize(sizei * 3, Nf);
343 }
344
345 for (label k = 0; k < Nf; k++)
346 {
347 List<Eigen::VectorXd> temp;
348 temp = field2EigenBC(fields[k]);
349
350 for (label i = 0; i < NBound; i++)
351 {
352 Out[i].col(k) = temp[i];
353 }
354 }
355
356 return Out;
357}
358
359template List<Eigen::MatrixXd> Foam2Eigen::PtrList2EigenBC(
360 PtrList<volVectorField>& fields, label Nfields);
361
362
363template <template <class> class PatchField, class GeoMesh>
364List<Eigen::MatrixXd> Foam2Eigen::PtrList2EigenBC(
365 PtrList<GeometricField<tensor, PatchField, GeoMesh>>&
366 fields,
367 label Nfields)
368{
369 label Nf;
370 M_Assert(Nfields <= fields.size(),
371 "The Number of requested fields cannot be bigger than the number of requested entries.");
372
373 if (Nfields == -1)
374 {
375 Nf = fields.size();
376 }
377 else
378 {
379 Nf = Nfields;
380 }
381
382 List<Eigen::MatrixXd> Out;
383 label NBound = fields[0].boundaryField().size();
384 Out.resize(NBound);
385
386 for (label i = 0; i < NBound; i++)
387 {
388 label sizei = fields[0].boundaryField()[i].size();
389 Out[i].resize(sizei * 9, Nf);
390 }
391
392 for (label k = 0; k < Nf; k++)
393 {
394 List<Eigen::VectorXd> temp;
395 temp = field2EigenBC(fields[k]);
396
397 for (label i = 0; i < NBound; i++)
398 {
399 Out[i].col(k) = temp[i];
400 }
401 }
402
403 return Out;
404}
405
406template List<Eigen::MatrixXd> Foam2Eigen::PtrList2EigenBC(
407 PtrList<volTensorField>& fields, label Nfields);
408
409template <template <class> class PatchField, class GeoMesh>
410GeometricField<tensor, PatchField, GeoMesh> Foam2Eigen::Eigen2field(
411 GeometricField<tensor, PatchField, GeoMesh>& field_in,
412 Eigen::VectorXd& eigen_vector, bool correctBC)
413{
414 GeometricField<tensor, PatchField, GeoMesh> field_out(field_in);
415
416 for (auto i = 0; i < field_out.size(); i++)
417 {
418 for (label j = 0; j < 9; j++)
419 {
420 field_out.ref()[i][j] = eigen_vector(i + field_out.size() * j);
421 }
422 }
423
424 if (correctBC)
425 {
426 field_out.correctBoundaryConditions();
427 }
428
429 return field_out;
430}
431
432template volTensorField Foam2Eigen::Eigen2field(
433 volTensorField& field_in, Eigen::VectorXd& eigen_vector, bool correctBC);
434
435template <template <class> class PatchField, class GeoMesh>
436GeometricField<vector, PatchField, GeoMesh> Foam2Eigen::Eigen2field(
437 GeometricField<vector, PatchField, GeoMesh>& field_in,
438 Eigen::VectorXd& eigen_vector, List<Eigen::VectorXd>& eigen_vector_boundary)
439{
440 GeometricField<vector, PatchField, GeoMesh> field_out(field_in);
441
442 for (auto i = 0; i < field_out.size(); i++)
443 {
444 field_out.ref()[i][0] = eigen_vector(i);
445 field_out.ref()[i][1] = eigen_vector(i + field_out.size());
446 field_out.ref()[i][2] = eigen_vector(i + field_out.size() * 2);
447 }
448
449 for(unsigned int id = 0; id < field_out.boundaryField().size(); id++)
450 {
451 unsigned int idBSize = field_out.boundaryField()[id].size();
452 for (unsigned int ith_bcell = 0; ith_bcell < idBSize; ith_bcell++)
453 {
454 ITHACAutilities::assignBC(field_out, id, eigen_vector_boundary[id]);
455 }
456 }
457
458 return field_out;
459}
460
461template volVectorField Foam2Eigen::Eigen2field(
462 volVectorField& field_in, Eigen::VectorXd& eigen_vector, List<Eigen::VectorXd>& eigen_vector_boundary);
463
464template<template<class> class PatchField, class GeoMesh>
465GeometricField<scalar, PatchField, GeoMesh> Foam2Eigen::Eigen2field(
466 GeometricField<scalar, PatchField, GeoMesh>& field_in,
467 Eigen::VectorXd& eigen_vector, List<Eigen::VectorXd>& eigen_vector_boundary)
468{
469 GeometricField<scalar, PatchField, GeoMesh> field_out(field_in);
470
471 for (auto i = 0; i < field_out.size(); i++)
472 {
473 field_out.ref()[i] = eigen_vector(i);
474 }
475
476 for(unsigned int id = 0; id < field_out.boundaryField().size(); id++)
477 {
478 for (unsigned int ith_bcell = 0; ith_bcell < field_out.boundaryField()[id].size(); ith_bcell++)
479 {
480 ITHACAutilities::assignBC(field_out, id, eigen_vector_boundary[id]);
481 }
482 }
483
484 return field_out;
485}
486
487template volScalarField Foam2Eigen::Eigen2field(
488 volScalarField& field_in, Eigen::VectorXd& eigen_vector, List<Eigen::VectorXd>& eigen_vector_boundary);
489
490template<template<class> class PatchField, class GeoMesh>
491GeometricField<vector, PatchField, GeoMesh> Foam2Eigen::Eigen2field(
492 GeometricField<vector, PatchField, GeoMesh>& field_in,
493 Eigen::VectorXd& eigen_vector, bool correctBC)
494{
495 GeometricField<vector, PatchField, GeoMesh> field_out(field_in);
496
497 for (auto i = 0; i < field_out.size(); i++)
498 {
499 for (label j = 0; j < 3; j++)
500 {
501 field_out.ref()[i][j] = eigen_vector(i + field_out.size() * j);
502 }
503 }
504
505 if (correctBC)
506 {
507 field_out.correctBoundaryConditions();
508 }
509
510 return field_out;
511}
512
513template volVectorField Foam2Eigen::Eigen2field(
514 volVectorField& field_in, Eigen::VectorXd& eigen_vector, bool correctBC);
515
516template <template <class> class PatchField, class GeoMesh>
517GeometricField<scalar, PatchField, GeoMesh> Foam2Eigen::Eigen2field(
518 GeometricField<scalar, PatchField, GeoMesh>& field_in,
519 Eigen::VectorXd& eigen_vector, bool correctBC)
520{
521 GeometricField<scalar, PatchField, GeoMesh> field_out(field_in);
522
523 for (auto i = 0; i < field_out.size(); i++)
524 {
525 field_out.ref()[i] = eigen_vector(i);
526 }
527
528 return field_out;
529}
530
531template surfaceScalarField Foam2Eigen::Eigen2field(
532 surfaceScalarField& field_in,
533 Eigen::VectorXd& eigen_vector,
534 bool correctBC);
535
536template <>
537volScalarField Foam2Eigen::Eigen2field(
538 volScalarField& field_in, Eigen::VectorXd& eigen_vector, bool correctBC)
539{
540 GeometricField<scalar, fvPatchField, volMesh> field_out(field_in);
541
542 for (auto i = 0; i < field_out.size(); i++)
543 {
544 field_out.ref()[i] = eigen_vector(i);
545 }
546
547 if (correctBC)
548 {
549 field_out.correctBoundaryConditions();
550 }
551
552 return field_out;
553}
554
555template <>
556Field<scalar> Foam2Eigen::Eigen2field(
557 Field<scalar>& field, Eigen::MatrixXd& matrix, bool correctBC)
558{
559 label sizeBC = field.size();
560 M_Assert(matrix.cols() == 1,
561 "The number of columns of the Input members is not correct, it should be 1");
562
563 if (matrix.rows() == 1)
564 {
565 Eigen::MatrixXd new_matrix = matrix.replicate(sizeBC, 1);
566 matrix.conservativeResize(sizeBC, 1);
567 matrix = new_matrix;
568 }
569
570 std::string message = "The input Eigen::MatrixXd has size " + name(
571 matrix.rows()) +
572 ". It should have the same size of the Field, i.e. " +
573 name(sizeBC);
574 M_Assert(matrix.rows() == sizeBC, message.c_str());
575
576 for (auto i = 0; i < sizeBC; i++)
577 {
578 field[i] = matrix(i, 0);
579 }
580
581 return field;
582}
583
584template <>
585Field<vector> Foam2Eigen::Eigen2field(
586 Field<vector>& field, Eigen::MatrixXd& matrix, bool correctBC)
587{
588 label sizeBC = field.size();
589 M_Assert(matrix.cols() == 3,
590 "The number of columns of the Input members is not correct, it should be 1");
591
592 if (matrix.rows() == 1)
593 {
594 Eigen::MatrixXd new_matrix = matrix.replicate(sizeBC, 1);
595 matrix.conservativeResize(sizeBC, 3);
596 matrix = new_matrix;
597 }
598
599 std::string message = "The input Eigen::MatrixXd has size " + name(
600 matrix.rows()) +
601 ". It should have the same size of the Field, i.e. " +
602 name(sizeBC);
603 M_Assert(matrix.rows() == sizeBC, message.c_str());
604
605 for (auto i = 0; i < sizeBC; i++)
606 {
607 for (label j = 0; j < 3; j++)
608 {
609 field[i][j] = matrix(i, j);
610 }
611 }
612
613 return field;
614}
615
616template <>
617Field<tensor> Foam2Eigen::Eigen2field(
618 Field<tensor>& field, Eigen::MatrixXd& matrix, bool correctBC)
619{
620 label sizeBC = field.size();
621 M_Assert(matrix.cols() == 9,
622 "The number of columns of the Input members is not correct, it should be 1");
623
624 if (matrix.rows() == 1)
625 {
626 Eigen::MatrixXd new_matrix = matrix.replicate(sizeBC, 1);
627 matrix.conservativeResize(sizeBC, 9);
628 matrix = new_matrix;
629 }
630
631 std::string message = "The input Eigen::MatrixXd has size " + name(
632 matrix.rows()) +
633 ". It should have the same size of the Field, i.e. " +
634 name(sizeBC);
635 M_Assert(matrix.rows() == sizeBC, message.c_str());
636
637 for (auto i = 0; i < sizeBC; i++)
638 {
639 for (label j = 0; j < 9; j++)
640 {
641 field[i][j] = matrix(i, j);
642 }
643 }
644
645 return field;
646}
647
648template <class Type, template <class> class PatchField, class GeoMesh>
649Eigen::MatrixXd Foam2Eigen::PtrList2Eigen(
650 PtrList<GeometricField<Type, PatchField, GeoMesh>>& fields,
651 label Nfields)
652{
653 label Nf;
654 M_Assert(Nfields <= fields.size(),
655 "The Number of requested fields cannot be bigger than the number of requested entries.");
656
657 if (Nfields == -1)
658 {
659 Nf = fields.size();
660 }
661 else
662 {
663 Nf = Nfields;
664 }
665
666 Eigen::MatrixXd out;
667 label nrows = (field2Eigen(fields[0])).rows();
668 out.resize(nrows, Nf);
669
670 for (label k = 0; k < Nf; k++)
671 {
672 out.col(k) = field2Eigen(fields[k]);
673 }
674
675 return out;
676}
677
678template Eigen::MatrixXd
679Foam2Eigen::PtrList2Eigen<scalar, fvPatchField, volMesh>(PtrList<volScalarField>&
680 fields,
681 label Nfields);
682template Eigen::MatrixXd Foam2Eigen::PtrList2Eigen(PtrList<surfaceScalarField>&
683 fields,
684 label Nfields);
685template Eigen::MatrixXd
686Foam2Eigen::PtrList2Eigen<vector, fvPatchField, volMesh>(PtrList<volVectorField>&
687 fields,
688 label Nfields);
689template Eigen::MatrixXd
690Foam2Eigen::PtrList2Eigen<tensor, fvPatchField, volMesh>(PtrList<volTensorField>&
691 fields,
692 label Nfields);
693
694template <>
695void Foam2Eigen::fvMatrix2Eigen(fvMatrix<scalar> foam_matrix,
696 Eigen::MatrixXd& A,
697 Eigen::VectorXd& b)
698{
699 label sizeA = foam_matrix.diag().size();
700 A.setZero(sizeA, sizeA);
701 b.setZero(sizeA);
702
703 for (auto i = 0; i < sizeA; i++)
704 {
705 A(i, i) = foam_matrix.diag()[i];
706 b(i, 0) = foam_matrix.source()[i];
707 }
708
709 const lduAddressing& addr = foam_matrix.lduAddr();
710 const labelList& lowerAddr = addr.lowerAddr();
711 const labelList& upperAddr = addr.upperAddr();
712 forAll(lowerAddr, i)
713 {
714 A(lowerAddr[i], upperAddr[i]) = foam_matrix.upper()[i];
715 A(upperAddr[i], lowerAddr[i]) = foam_matrix.lower()[i];
716 }
717 forAll(foam_matrix.psi().boundaryField(), I)
718 {
719 const fvPatch& ptch = foam_matrix.psi().boundaryField()[I].patch();
720 forAll(ptch, J)
721 {
722 label w = ptch.faceCells()[J];
723 const double intern = foam_matrix.internalCoeffs()[I][J];
724 A(w, w) += intern;
725 b(w, 0) += foam_matrix.boundaryCoeffs()[I][J];
726 }
727 }
728}
729
730template <>
731void Foam2Eigen::fvMatrix2Eigen(fvMatrix<vector> foam_matrix,
732 Eigen::MatrixXd& A,
733 Eigen::VectorXd& b)
734{
735 label sizeA = foam_matrix.diag().size();
736 A.resize(sizeA * 3, sizeA * 3);
737 b.resize(sizeA * 3);
738
739 for (auto i = 0; i < sizeA; i++)
740 {
741 A(i, i) = foam_matrix.diag()[i];
742 A(sizeA + i, sizeA + i) = foam_matrix.diag()[i];
743 A(2 * sizeA + i, 2 * sizeA + i) = foam_matrix.diag()[i];
744 b(i) = foam_matrix.source()[i][0];
745 b(sizeA + i) = foam_matrix.source()[i][1];
746 b(2 * sizeA + i) = foam_matrix.source()[i][2];
747 }
748
749 const lduAddressing& addr = foam_matrix.lduAddr();
750 const labelList& lowerAddr = addr.lowerAddr();
751 const labelList& upperAddr = addr.upperAddr();
752 forAll(lowerAddr, i)
753 {
754 A(lowerAddr[i], upperAddr[i]) = foam_matrix.upper()[i];
755 A(lowerAddr[i] + sizeA, upperAddr[i] + sizeA) = foam_matrix.upper()[i];
756 A(lowerAddr[i] + sizeA * 2, upperAddr[i] + sizeA * 2) = foam_matrix.upper()[i];
757 A(upperAddr[i], lowerAddr[i]) = foam_matrix.lower()[i];
758 A(upperAddr[i] + sizeA, lowerAddr[i] + sizeA) = foam_matrix.lower()[i];
759 A(upperAddr[i] + sizeA * 2, lowerAddr[i] + sizeA * 2) = foam_matrix.lower()[i];
760 }
761 forAll(foam_matrix.psi().boundaryField(), I)
762 {
763 const fvPatch& ptch = foam_matrix.psi().boundaryField()[I].patch();
764 forAll(ptch, J)
765 {
766 label w = ptch.faceCells()[J];
767 A(w, w) += foam_matrix.internalCoeffs()[I][J][0];
768 A(w + sizeA, w + sizeA) += foam_matrix.internalCoeffs()[I][J][1];
769 A(w + sizeA * 2, w + sizeA * 2) += foam_matrix.internalCoeffs()[I][J][2];
770 b(w) += foam_matrix.boundaryCoeffs()[I][J][0];
771 b(w + sizeA) += foam_matrix.boundaryCoeffs()[I][J][1];
772 b(w + sizeA * 2) += foam_matrix.boundaryCoeffs()[I][J][2];
773 }
774 }
775}
776
777template <>
778void Foam2Eigen::fvMatrix2Eigen(fvMatrix<scalar> foam_matrix,
779 Eigen::SparseMatrix<double>& A, Eigen::VectorXd& b)
780{
781 label sizeA = foam_matrix.diag().size();
782 label nel = foam_matrix.diag().size() + foam_matrix.upper().size() +
783 foam_matrix.lower().size();
784 A.resize(sizeA, sizeA);
785 b.resize(sizeA);
786 A.reserve(nel);
787 typedef Eigen::Triplet<double> Trip;
788 std::vector<Trip> tripletList;
789 tripletList.reserve(nel);
790
791 for (auto i = 0; i < sizeA; i++)
792 {
793 tripletList.push_back(Trip(i, i, foam_matrix.diag()[i]));
794 b(i, 0) = foam_matrix.source()[i];
795 }
796
797 const lduAddressing& addr = foam_matrix.lduAddr();
798 const labelList& lowerAddr = addr.lowerAddr();
799 const labelList& upperAddr = addr.upperAddr();
800 forAll(lowerAddr, i)
801 {
802 tripletList.push_back(Trip(lowerAddr[i], upperAddr[i], foam_matrix.upper()[i]));
803 tripletList.push_back(Trip(upperAddr[i], lowerAddr[i], foam_matrix.lower()[i]));
804 }
805 forAll(foam_matrix.psi().boundaryField(), I)
806 {
807 const fvPatch& ptch = foam_matrix.psi().boundaryField()[I].patch();
808 forAll(ptch, J)
809 {
810 label w = ptch.faceCells()[J];
811 tripletList.push_back(Trip(w, w, foam_matrix.internalCoeffs()[I][J]));
812 b(w, 0) += foam_matrix.boundaryCoeffs()[I][J];
813 }
814 }
815 A.setFromTriplets(tripletList.begin(), tripletList.end());
816}
817
818template <>
819void Foam2Eigen::fvMatrix2Eigen(fvMatrix<vector> foam_matrix,
820 Eigen::SparseMatrix<double>& A, Eigen::VectorXd& b)
821{
822 label sizeA = foam_matrix.diag().size();
823 label nel = foam_matrix.diag().size() + foam_matrix.upper().size() +
824 foam_matrix.lower().size();
825 A.resize(sizeA * 3, sizeA * 3);
826 A.reserve(nel * 3);
827 b.resize(sizeA * 3);
828 typedef Eigen::Triplet<double> Trip;
829 std::vector<Trip> tripletList;
830 tripletList.reserve(nel * 3);
831
832 for (auto i = 0; i < sizeA; i++)
833 {
834 tripletList.push_back(Trip(i, i, foam_matrix.diag()[i]));
835 tripletList.push_back(Trip(sizeA + i, sizeA + i, foam_matrix.diag()[i]));
836 tripletList.push_back(Trip(2 * sizeA + i, 2 * sizeA + i,
837 foam_matrix.diag()[i]));
838 b(i) = foam_matrix.source()[i][0];
839 b(sizeA + i) = foam_matrix.source()[i][1];
840 b(2 * sizeA + i) = foam_matrix.source()[i][2];
841 }
842
843 const lduAddressing& addr = foam_matrix.lduAddr();
844 const labelList& lowerAddr = addr.lowerAddr();
845 const labelList& upperAddr = addr.upperAddr();
846 forAll(lowerAddr, i)
847 {
848 tripletList.push_back(Trip(lowerAddr[i], upperAddr[i], foam_matrix.upper()[i]));
849 tripletList.push_back(Trip(lowerAddr[i] + sizeA, upperAddr[i] + sizeA,
850 foam_matrix.upper()[i]));
851 tripletList.push_back(Trip(lowerAddr[i] + sizeA * 2, upperAddr[i] + sizeA * 2,
852 foam_matrix.upper()[i]));
853 tripletList.push_back(Trip(upperAddr[i], lowerAddr[i], foam_matrix.lower()[i]));
854 tripletList.push_back(Trip(upperAddr[i] + sizeA, lowerAddr[i] + sizeA,
855 foam_matrix.lower()[i]));
856 tripletList.push_back(Trip(upperAddr[i] + sizeA * 2, lowerAddr[i] + sizeA * 2,
857 foam_matrix.lower()[i]));
858 }
859 forAll(foam_matrix.psi().boundaryField(), I)
860 {
861 const fvPatch& ptch = foam_matrix.psi().boundaryField()[I].patch();
862 forAll(ptch, J)
863 {
864 label w = ptch.faceCells()[J];
865 tripletList.push_back(Trip(w, w, foam_matrix.internalCoeffs()[I][J][0]));
866 tripletList.push_back(Trip(w + sizeA, w + sizeA,
867 foam_matrix.internalCoeffs()[I][J][1]));
868 tripletList.push_back(Trip(w + sizeA * 2, w + sizeA * 2,
869 foam_matrix.internalCoeffs()[I][J][2]));
870 b(w) += foam_matrix.boundaryCoeffs()[I][J][0];
871 b(w + sizeA) += foam_matrix.boundaryCoeffs()[I][J][1];
872 b(w + sizeA * 2) += foam_matrix.boundaryCoeffs()[I][J][2];
873 }
874 }
875 A.setFromTriplets(tripletList.begin(), tripletList.end());
876}
877
878template <>
879void Foam2Eigen::fvMatrix2EigenM(fvMatrix<scalar>& foam_matrix,
880 Eigen::MatrixXd& A)
881{
882 label sizeA = foam_matrix.diag().size();
883 A.setZero(sizeA, sizeA);
884
885 for (auto i = 0; i < sizeA; i++)
886 {
887 A(i, i) = foam_matrix.diag()[i];
888 }
889
890 const lduAddressing& addr = foam_matrix.lduAddr();
891 const labelList& lowerAddr = addr.lowerAddr();
892 const labelList& upperAddr = addr.upperAddr();
893 forAll(lowerAddr, i)
894 {
895 A(lowerAddr[i], upperAddr[i]) = foam_matrix.upper()[i];
896 A(upperAddr[i], lowerAddr[i]) = foam_matrix.lower()[i];
897 }
898 forAll(foam_matrix.psi().boundaryField(), I)
899 {
900 const fvPatch& ptch = foam_matrix.psi().boundaryField()[I].patch();
901 forAll(ptch, J)
902 {
903 label w = ptch.faceCells()[J];
904 A(w, w) += foam_matrix.internalCoeffs()[I][J];
905 }
906 }
907}
908
909template <>
910void Foam2Eigen::fvMatrix2EigenM(fvMatrix<scalar>& foam_matrix,
911 Eigen::SparseMatrix<double>& A)
912{
913 label sizeA = foam_matrix.diag().size();
914 label nel = foam_matrix.diag().size() + foam_matrix.upper().size() +
915 foam_matrix.lower().size();
916 A.resize(sizeA, sizeA);
917 A.reserve(nel);
918 typedef Eigen::Triplet<double> Trip;
919 std::vector<Trip> tripletList;
920 tripletList.reserve(nel);
921
922 for (auto i = 0; i < sizeA; i++)
923 {
924 tripletList.push_back(Trip(i, i, foam_matrix.diag()[i]));
925 }
926
927 const lduAddressing& addr = foam_matrix.lduAddr();
928 const labelList& lowerAddr = addr.lowerAddr();
929 const labelList& upperAddr = addr.upperAddr();
930 forAll(lowerAddr, i)
931 {
932 tripletList.push_back(Trip(lowerAddr[i], upperAddr[i], foam_matrix.upper()[i]));
933 tripletList.push_back(Trip(upperAddr[i], lowerAddr[i], foam_matrix.lower()[i]));
934 }
935 forAll(foam_matrix.psi().boundaryField(), I)
936 {
937 const fvPatch& ptch = foam_matrix.psi().boundaryField()[I].patch();
938 forAll(ptch, J)
939 {
940 label w = ptch.faceCells()[J];
941 tripletList.push_back(Trip(w, w, foam_matrix.internalCoeffs()[I][J]));
942 }
943 }
944 A.setFromTriplets(tripletList.begin(), tripletList.end());
945}
946
947template <>
948void Foam2Eigen::fvMatrix2EigenM(fvMatrix<vector>& foam_matrix,
949 Eigen::MatrixXd& A)
950{
951 label sizeA = foam_matrix.diag().size();
952 A.resize(sizeA * 3, sizeA * 3);
953
954 for (auto i = 0; i < sizeA; i++)
955 {
956 A(i, i) = foam_matrix.diag()[i];
957 A(sizeA + i, sizeA + i) = foam_matrix.diag()[i];
958 A(2 * sizeA + i, 2 * sizeA + i) = foam_matrix.diag()[i];
959 }
960
961 const lduAddressing& addr = foam_matrix.lduAddr();
962 const labelList& lowerAddr = addr.lowerAddr();
963 const labelList& upperAddr = addr.upperAddr();
964 forAll(lowerAddr, i)
965 {
966 A(lowerAddr[i], upperAddr[i]) = foam_matrix.upper()[i];
967 A(lowerAddr[i] + sizeA, upperAddr[i] + sizeA) = foam_matrix.upper()[i];
968 A(lowerAddr[i] + sizeA * 2, upperAddr[i] + sizeA * 2) = foam_matrix.upper()[i];
969 A(upperAddr[i], lowerAddr[i]) = foam_matrix.lower()[i];
970 A(upperAddr[i] + sizeA, lowerAddr[i] + sizeA) = foam_matrix.lower()[i];
971 A(upperAddr[i] + sizeA * 2, lowerAddr[i] + sizeA * 2) = foam_matrix.lower()[i];
972 }
973 forAll(foam_matrix.psi().boundaryField(), I)
974 {
975 const fvPatch& ptch = foam_matrix.psi().boundaryField()[I].patch();
976 forAll(ptch, J)
977 {
978 label w = ptch.faceCells()[J];
979 A(w, w) += foam_matrix.internalCoeffs()[I][J][0];
980 A(w + sizeA, w + sizeA) += foam_matrix.internalCoeffs()[I][J][1];
981 A(w + sizeA * 2, w + sizeA * 2) += foam_matrix.internalCoeffs()[I][J][2];
982 }
983 }
984}
985
986
987template <>
988void Foam2Eigen::fvMatrix2EigenM(fvMatrix<vector>& foam_matrix,
989 Eigen::SparseMatrix<double>& A)
990{
991 label sizeA = foam_matrix.diag().size();
992 label nel = foam_matrix.diag().size() + foam_matrix.upper().size() +
993 foam_matrix.lower().size();
994 A.resize(sizeA * 3, sizeA * 3);
995 A.reserve(nel * 3);
996 typedef Eigen::Triplet<double> Trip;
997 std::vector<Trip> tripletList;
998 tripletList.reserve(nel * 3);
999
1000 for (auto i = 0; i < sizeA; i++)
1001 {
1002 tripletList.push_back(Trip(i, i, foam_matrix.diag()[i]));
1003 tripletList.push_back(Trip(sizeA + i, sizeA + i, foam_matrix.diag()[i]));
1004 tripletList.push_back(Trip(2 * sizeA + i, 2 * sizeA + i,
1005 foam_matrix.diag()[i]));
1006 }
1007
1008 const lduAddressing& addr = foam_matrix.lduAddr();
1009 const labelList& lowerAddr = addr.lowerAddr();
1010 const labelList& upperAddr = addr.upperAddr();
1011 forAll(lowerAddr, i)
1012 {
1013 tripletList.push_back(Trip(lowerAddr[i], upperAddr[i], foam_matrix.upper()[i]));
1014 tripletList.push_back(Trip(lowerAddr[i] + sizeA, upperAddr[i] + sizeA,
1015 foam_matrix.upper()[i]));
1016 tripletList.push_back(Trip(lowerAddr[i] + sizeA * 2, upperAddr[i] + sizeA * 2,
1017 foam_matrix.upper()[i]));
1018 tripletList.push_back(Trip(upperAddr[i], lowerAddr[i], foam_matrix.lower()[i]));
1019 tripletList.push_back(Trip(upperAddr[i] + sizeA, lowerAddr[i] + sizeA,
1020 foam_matrix.lower()[i]));
1021 tripletList.push_back(Trip(upperAddr[i] + sizeA * 2, lowerAddr[i] + sizeA * 2,
1022 foam_matrix.lower()[i]));
1023 }
1024 forAll(foam_matrix.psi().boundaryField(), I)
1025 {
1026 const fvPatch& ptch = foam_matrix.psi().boundaryField()[I].patch();
1027 forAll(ptch, J)
1028 {
1029 label w = ptch.faceCells()[J];
1030 tripletList.push_back(Trip(w, w, foam_matrix.internalCoeffs()[I][J][0]));
1031 tripletList.push_back(Trip(w + sizeA, w + sizeA,
1032 foam_matrix.internalCoeffs()[I][J][1]));
1033 tripletList.push_back(Trip(w + sizeA * 2, w + sizeA * 2,
1034 foam_matrix.internalCoeffs()[I][J][2]));
1035 }
1036 }
1037 A.setFromTriplets(tripletList.begin(), tripletList.end());
1038}
1039
1040template <>
1041void Foam2Eigen::fvMatrix2EigenV(fvMatrix<scalar>& foam_matrix,
1042 Eigen::VectorXd& b)
1043{
1044 label sizeA = foam_matrix.diag().size();
1045 b.setZero(sizeA);
1046
1047 for (auto i = 0; i < sizeA; i++)
1048 {
1049 b(i, 0) = foam_matrix.source()[i];
1050 }
1051
1052 forAll(foam_matrix.psi().boundaryField(), I)
1053 {
1054 const fvPatch& ptch = foam_matrix.psi().boundaryField()[I].patch();
1055 forAll(ptch, J)
1056 {
1057 label w = ptch.faceCells()[J];
1058 b(w, 0) += foam_matrix.boundaryCoeffs()[I][J];
1059 }
1060 }
1061}
1062
1063template <>
1064void Foam2Eigen::fvMatrix2EigenV(fvMatrix<vector>& foam_matrix,
1065 Eigen::VectorXd& b)
1066{
1067 label sizeA = foam_matrix.diag().size();
1068 b.resize(sizeA * 3);
1069
1070 for (auto i = 0; i < sizeA; i++)
1071 {
1072 b(i) = foam_matrix.source()[i][0];
1073 b(sizeA + i) = foam_matrix.source()[i][1];
1074 b(2 * sizeA + i) = foam_matrix.source()[i][2];
1075 }
1076
1077 forAll(foam_matrix.psi().boundaryField(), I)
1078 {
1079 const fvPatch& ptch = foam_matrix.psi().boundaryField()[I].patch();
1080 forAll(ptch, J)
1081 {
1082 label w = ptch.faceCells()[J];
1083 b(w) += foam_matrix.boundaryCoeffs()[I][J][0];
1084 b(w + sizeA) += foam_matrix.boundaryCoeffs()[I][J][1];
1085 b(w + sizeA * 2) += foam_matrix.boundaryCoeffs()[I][J][2];
1086 }
1087 }
1088}
1089
1090template <class Type, template <class> class PatchField, class GeoMesh>
1091Eigen::VectorXd Foam2Eigen::projectField(
1092 GeometricField<Type, PatchField, GeoMesh>& field,
1093 PtrList<GeometricField<Type, PatchField, GeoMesh>>& modes,
1094 label Nmodes)
1095{
1096 Eigen::VectorXd fr;
1097 Eigen::MatrixXd Eig_Modes = PtrList2Eigen(modes, Nmodes);
1098 Eigen::VectorXd f = Foam2Eigen::field2Eigen(field);
1099 Eigen::VectorXd Volumes = field2Eigen(modes[0].mesh());
1100 Eigen::MatrixXd VolumesN(Volumes.rows(), 1);
1101 VolumesN = Volumes;
1102
1103 if (Volumes.rows() != Eig_Modes.rows())
1104 {
1105 VolumesN.resize(Eig_Modes.rows(), 1);
1106 VolumesN.col(0).segment(0, Volumes.rows()) = Volumes;
1107 VolumesN.col(0).segment(Volumes.rows() + 1, Volumes.rows() * 2) = Volumes;
1108 VolumesN.col(0).segment(Volumes.rows() * 2 + 1, Volumes.rows() * 3) = Volumes;
1109 }
1110
1111 fr = Eig_Modes.transpose() * (f.cwiseProduct(VolumesN));
1112 return fr;
1113}
1114
1115template <class Type, template <class> class PatchField, class GeoMesh>
1116std::tuple<Eigen::MatrixXd, Eigen::VectorXd> Foam2Eigen::projectFvMatrix(
1117 fvMatrix<Type>& matrix,
1118 PtrList<GeometricField<Type, PatchField, GeoMesh>>& modes, label Nmodes)
1119{
1120 Eigen::SparseMatrix<double> A;
1121 Eigen::MatrixXd Ar;
1122 Eigen::VectorXd b;
1123 Eigen::VectorXd br;
1124 Eigen::MatrixXd Eig_Modes = PtrList2Eigen(modes, Nmodes);
1125 Foam2Eigen::fvMatrix2Eigen(matrix, A, b);
1126 Eigen::VectorXd Volumes = field2Eigen(modes[0].mesh());
1127 Eigen::MatrixXd VolumesN(Volumes.rows(), Nmodes);
1128
1129 if (Volumes.rows() != Eig_Modes.rows())
1130 {
1131 VolumesN.resize(Eig_Modes.rows(), Nmodes);
1132 }
1133
1134 if (Volumes.rows() == Eig_Modes.rows())
1135 {
1136 for (auto i = 0; i < Nmodes; i++)
1137 {
1138 VolumesN.col(i) = Volumes;
1139 }
1140 }
1141 else
1142 {
1143 for (auto i = 0; i < Nmodes; i++)
1144 {
1145 VolumesN.col(i).segment(0, Volumes.rows()) = Volumes;
1146 VolumesN.col(i).segment(Volumes.rows() + 1, Volumes.rows() * 2) = Volumes;
1147 VolumesN.col(i).segment(Volumes.rows() * 2 + 1, Volumes.rows() * 3) = Volumes;
1148 }
1149 }
1150
1151 Ar = Eig_Modes.transpose() * A * Eig_Modes;
1152 br = Eig_Modes.transpose() * b;
1153 std::tuple<Eigen::MatrixXd, Eigen::VectorXd> tupla;
1154 tupla = std::make_tuple(Ar, br);
1155 return tupla;
1156}
1157
1158template <class Type, template <class> class PatchField, class GeoMesh>
1159Eigen::MatrixXd Foam2Eigen::MassMatrix(
1160 PtrList<GeometricField<Type, PatchField, GeoMesh>>& modes, label Nmodes)
1161{
1162 Eigen::MatrixXd Mr;
1163 Eigen::MatrixXd Eig_Modes = PtrList2Eigen(modes, Nmodes);
1164 Eigen::VectorXd Volumes = field2Eigen(modes[0].mesh());
1165 Eigen::MatrixXd VolumesN(Volumes.rows(), Nmodes);
1166
1167 if (Volumes.rows() != Eig_Modes.rows())
1168 {
1169 VolumesN.resize(Eig_Modes.rows(), Nmodes);
1170 }
1171
1172 if (Volumes.rows() == Eig_Modes.rows())
1173 {
1174 for (auto i = 0; i < Nmodes; i++)
1175 {
1176 VolumesN.col(i) = Volumes;
1177 }
1178 }
1179 else
1180 {
1181 for (auto i = 0; i < Nmodes; i++)
1182 {
1183 VolumesN.col(i).segment(0, Volumes.rows()) = Volumes;
1184 VolumesN.col(i).segment(Volumes.rows() + 1, Volumes.rows() * 2) = Volumes;
1185 VolumesN.col(i).segment(Volumes.rows() * 2 + 1, Volumes.rows() * 3) = Volumes;
1186 }
1187 }
1188
1189 Mr = Eig_Modes.transpose() * (Eig_Modes.cwiseProduct(VolumesN));
1190 return Mr;
1191}
1192
1193template <class Type>
1194std::tuple<List<Eigen::SparseMatrix<double>>, List<Eigen::VectorXd>>
1195 Foam2Eigen::LFvMatrix2LSM(PtrList<fvMatrix<Type>>& MatrixList)
1196{
1197 List<Eigen::SparseMatrix<double>> SM_list;
1198 List<Eigen::VectorXd> V_list;
1199 label LSize = MatrixList.size();
1200 SM_list.resize(LSize);
1201 V_list.resize(LSize);
1202 Eigen::SparseMatrix<double> A;
1203 Eigen::VectorXd b;
1204
1205 for (label j = 0; j < LSize; j++)
1206 {
1207 fvMatrix2Eigen(MatrixList[j], A, b);
1208 SM_list[j] = A;
1209 V_list[j] = b;
1210 }
1211
1212 std::tuple<List<Eigen::SparseMatrix<double>>, List<Eigen::VectorXd>> tupla;
1213 tupla = std::make_tuple(SM_list, V_list);
1214 return tupla;
1215}
1216
1217template std::tuple<List<Eigen::SparseMatrix<double>>, List<Eigen::VectorXd>>
1218Foam2Eigen::LFvMatrix2LSM(PtrList<fvMatrix<scalar>>& MatrixList);
1219template std::tuple<List<Eigen::SparseMatrix<double>>, List<Eigen::VectorXd>>
1220Foam2Eigen::LFvMatrix2LSM(PtrList<fvMatrix<vector>>& MatrixList);
1221
1222template <class type_matrix>
1223Eigen::Matrix<type_matrix, Eigen::Dynamic, Eigen::Dynamic>
1224Foam2Eigen::List2EigenMatrix(List<type_matrix> list)
1225{
1226 Eigen::Matrix<type_matrix, Eigen::Dynamic, Eigen::Dynamic> matrix(list.size(),
1227 1);
1228
1229 for (label i = 0; i < matrix.rows(); i++)
1230 {
1231 matrix(i, 0) = list[i];
1232 }
1233
1234 return matrix;
1235}
1236
1237template Eigen::Matrix<int, Eigen::Dynamic, Eigen::Dynamic>
1238Foam2Eigen::List2EigenMatrix(List<int> list);
1239template Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>
1240Foam2Eigen::List2EigenMatrix(List<double> list);
1241
1242template <class type_matrix>
1243List<type_matrix> Foam2Eigen::EigenMatrix2List(
1244 Eigen::Matrix<type_matrix, Eigen::Dynamic, Eigen::Dynamic> matrix)
1245{
1246 if (matrix.cols() == 1)
1247 {
1248 List<type_matrix> list(matrix.rows());
1249
1250 for (label i = 0; i < matrix.rows(); i++)
1251 {
1252 list[i] = matrix(i, 0);
1253 }
1254
1255 return list;
1256 }
1257 else
1258 {
1259 Info << "Foam2Eigen::EigenMatrix2List only accepts matrices with 1 column, exiting"
1260 << endl;
1261 exit(11);
1262 }
1263}
1264
1265template List<int> Foam2Eigen::EigenMatrix2List(
1266 Eigen::Matrix<int, Eigen::Dynamic, Eigen::Dynamic> matrix);
1267template List<double> Foam2Eigen::EigenMatrix2List(
1268 Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> matrix);
1269
1270template <>
1271Eigen::MatrixXd Foam2Eigen::field2Eigen(const List<vector>& field)
1272{
1273 Eigen::MatrixXd out;
1274 out.resize(label(field.size() * 3), 1);
1275
1276 for (label l = 0; l < field.size(); l++)
1277 {
1278 for (label j = 0; j < 3; j++)
1279 {
1280 out(j * field.size() + l, 0) = field[l][j];
1281 }
1282 }
1283
1284 return out;
1285}
1286
1287template <>
1288Eigen::MatrixXd Foam2Eigen::field2Eigen(const List<scalar>& field)
1289{
1290 Eigen::MatrixXd out;
1291 out.resize(label(field.size()), 1);
1292
1293 for (label l = 0; l < field.size(); l++)
1294 {
1295 out(l, 0) = field[l];
1296 }
1297
1298 return out;
1299}
forAll
forAll(example_CG.gList, solutionI)
Definition CGtest.H:21
Foam2Eigen.H
Header file of the Foam2Eigen class.
mesh
Foam::fvMesh & mesh
Definition createMesh.H:47
M_Assert
#define M_Assert(Expr, Msg)
Definition ITHACAassert.H:46
Foam2Eigen::Eigen2field
static GeometricField< scalar, PatchField, GeoMesh > Eigen2field(GeometricField< scalar, PatchField, GeoMesh > &field, Eigen::VectorXd &eigen_vector, bool correctBC=true)
Convert a vector in Eigen format into an OpenFOAM scalar GeometricField.
Definition Foam2Eigen.C:517
Foam2Eigen::field2Eigen
static Eigen::VectorXd field2Eigen(GeometricField< tensor, PatchField, GeoMesh > &field)
Convert a vector OpenFOAM field into an Eigen Vector.
Definition Foam2Eigen.C:40
Foam2Eigen::field2EigenBC
static List< Eigen::VectorXd > field2EigenBC(GeometricField< scalar, PatchField, GeoMesh > &field)
Convert an OpenFOAM scalar field to a List of Eigen Vectors, one for each boundary.
Definition Foam2Eigen.C:244
Foam2Eigen::field2EigenMap
static Eigen::Map< Eigen::MatrixXd > field2EigenMap(GeometricField< scalar, PatchField, GeoMesh > &field)
Convert a scalar OpenFOAM field into an Eigen Vector.
Definition Foam2Eigen.C:97
Foam2Eigen::LFvMatrix2LSM
static std::tuple< List< Eigen::SparseMatrix< double > >, List< Eigen::VectorXd > > LFvMatrix2LSM(PtrList< fvMatrix< Type > > &MatrixList)
Convert a PtrList of OpenFOAM fvMatrix into a tuple of lists of Eigen Sparse Matrices and source vect...
Definition Foam2Eigen.C:1195
Foam2Eigen::List2EigenMatrix
static Eigen::Matrix< type_matrix, Eigen::Dynamic, Eigen::Dynamic > List2EigenMatrix(List< type_matrix > list)
Convert a Foam List into an Eigen matrix with one column.
Definition Foam2Eigen.C:1224
Foam2Eigen::projectField
static Eigen::VectorXd projectField(GeometricField< Type, PatchField, GeoMesh > &field, PtrList< GeometricField< Type, PatchField, GeoMesh > > &modes, label Nmodes)
Perform the projection of an OpenFOAM field onto a set of modes using the Eigen matrix multiplication...
Definition Foam2Eigen.C:1091
Foam2Eigen::MassMatrix
static Eigen::MatrixXd MassMatrix(PtrList< GeometricField< Type, PatchField, GeoMesh > > &modes, label Nmodes)
Obtain the Mass Matrix from a list of reduced basis.
Definition Foam2Eigen.C:1159
Foam2Eigen::EigenMatrix2List
static List< type_matrix > EigenMatrix2List(Eigen::Matrix< type_matrix, Eigen::Dynamic, Eigen::Dynamic > matrix)
Convert an Eigen matrix with one column into a Foam List.
Definition Foam2Eigen.C:1243
Foam2Eigen::fvMatrix2EigenV
static void fvMatrix2EigenV(fvMatrix< type_foam_matrix > &foam_matrix, type_B &b)
Convert a ldu OpenFOAM matrix into a source vector b.
Foam2Eigen::PtrList2EigenBC
static List< Eigen::MatrixXd > PtrList2EigenBC(PtrList< GeometricField< scalar, PatchField, GeoMesh > > &fields, label Nfields=-1)
Convert an OpenFOAM scalar field to a List of Eigen Vectors, one for each boundary.
Definition Foam2Eigen.C:268
Foam2Eigen::PtrList2Eigen
static Eigen::MatrixXd PtrList2Eigen(PtrList< GeometricField< Type, PatchField, GeoMesh > > &fields, label Nfields=-1)
Convert a PtrList of snapshots to Eigen matrix (only internal field)
Definition Foam2Eigen.C:649
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.
Foam2Eigen::projectFvMatrix
static std::tuple< Eigen::MatrixXd, Eigen::VectorXd > projectFvMatrix(fvMatrix< Type > &matrix, PtrList< GeometricField< Type, PatchField, GeoMesh > > &modes, label Nmodes)
Obtain the Mass Matrix from a list of reduced basis.
Definition Foam2Eigen.C:1116
Foam2Eigen::field2EigenMapBC
static Eigen::Map< Eigen::MatrixXd > field2EigenMapBC(GeometricField< scalar, PatchField, GeoMesh > &field, int BC_index)
Map a scalar OpenFOAM field boundary into an Eigen Matrix.
Definition Foam2Eigen.C:105
Foam2Eigen::fvMatrix2EigenM
static void fvMatrix2EigenM(fvMatrix< type_foam_matrix > &foam_matrix, type_A &A)
Convert a ldu OpenFOAM matrix into a Eigen Matrix A.
A
volScalarField & A
Definition createConstants.H:144
ITHACAutilities::assignBC
void assignBC(GeometricField< scalar, fvPatchField, volMesh > &s, label BC_ind, double value)
Assign uniform Boundary Condition to a volScalarField.
Definition ITHACAassign.C:143
i
label i
Definition pEqn.H:46
Generated by   doxygen 1.11.0