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