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