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