From f12399c39b8cdc2ed5e6ad8ee6d535af88710366 Mon Sep 17 00:00:00 2001
From: Kutalmis Bercin <kutalmis.bercin@esi-group.com>
Date: Mon, 28 Feb 2022 13:50:05 +0000
Subject: [PATCH] ENH: DMD: add multi-patch input functionality

---
 src/functionObjects/field/DMD/DMD.C           | 24 +++++---
 src/functionObjects/field/DMD/DMD.H           | 22 ++++---
 .../field/DMD/DMDModels/derived/STDMD/STDMD.C |  9 ++-
 .../field/DMD/DMDModels/derived/STDMD/STDMD.H | 12 ++--
 .../DMDModels/derived/STDMD/STDMDTemplates.C  | 61 +++++++++++--------
 src/functionObjects/field/DMD/DMDTemplates.C  | 35 ++++++-----
 6 files changed, 101 insertions(+), 62 deletions(-)

diff --git a/src/functionObjects/field/DMD/DMD.C b/src/functionObjects/field/DMD/DMD.C
index d419543586..e67ba80d61 100644
--- a/src/functionObjects/field/DMD/DMD.C
+++ b/src/functionObjects/field/DMD/DMD.C
@@ -91,22 +91,21 @@ void Foam::functionObjects::DMD::initialise()
 {
     const label nComps = nComponents(fieldName_);
 
-    if (patch_.empty())
+    if (patches_.empty())
     {
         nSnap_ = nComps*mesh_.nCells();
     }
     else
     {
-        const label patchi = mesh_.boundaryMesh().findPatchID(patch_);
+        const labelList patchis
+        (
+            mesh_.boundaryMesh().patchSet(patches_).sortedToc()
+        );
 
-        if (patchi < 0)
+        for (const label patchi : patchis)
         {
-            FatalErrorInFunction
-                << "Cannot find patch " << patch_
-                << exit(FatalError);
+            nSnap_ += nComps*(mesh_.C().boundaryField()[patchi]).size();
         }
-
-        nSnap_ = nComps*(mesh_.C().boundaryField()[patchi]).size();
     }
 
     const label nSnapTotal = returnReduce(nSnap_, sumOp<label>());
@@ -141,8 +140,15 @@ Foam::functionObjects::DMD::DMD
     fvMeshFunctionObject(name, runTime, dict),
     DMDModelPtr_(DMDModel::New(mesh_, name, dict)),
     z_(),
+    patches_
+    (
+        dict.getOrDefault<wordRes>
+        (
+            "patches",
+            dict.found("patch") ? wordRes(1,dict.get<word>("patch")) : wordRes()
+        )
+    ),
     fieldName_(dict.get<word>("field")),
-    patch_(dict.getOrDefault<word>("patch", word::null)),
     nSnap_(0),
     step_(0)
 {
diff --git a/src/functionObjects/field/DMD/DMD.H b/src/functionObjects/field/DMD/DMD.H
index 669e7f2e5c..acabb983d0 100644
--- a/src/functionObjects/field/DMD/DMD.H
+++ b/src/functionObjects/field/DMD/DMD.H
@@ -87,16 +87,21 @@ Usage
     \verbatim
     DMD1
     {
-        // Mandatory entries (unmodifiable)
+        // Mandatory entries
         type                DMD;
         libs                (fieldFunctionObjects);
         DMDModel            <DMDModel>;
-        field               <inputField>;
+        field               <word>;
 
-        // Optional entries (unmodifiable)
-        patch               <patchName>;
+        // Optional entries
 
-        // Mandatory/Optional (inherited) entries
+            // Option-1
+            patch               <word>;
+
+            // Option-2
+            patches             (<wordRes>);
+
+        // Inherited entries
         ...
     }
     \endverbatim
@@ -109,6 +114,7 @@ Usage
       DMDModel     | Name of specified DMD model        | word | yes  | -
       field        | Name of operand field              | word | yes  | -
       patch        | Name of operand patch              | word | no   | null
+      patches      | Names of operand patches        | wordRes | no   | null
     \endtable
 
     Options for the \c DMDModel entry:
@@ -181,12 +187,12 @@ class DMD
         //  time step, an augmented snapshot is constructed of two snapshots
         RMatrix z_;
 
+        //- Names of operand patches
+        const wordRes patches_;
+
         //- Name of operand field
         const word fieldName_;
 
-        //- Name of operand patch
-        const word patch_;
-
         //- Number of elements in a snapshot
         label nSnap_;
 
diff --git a/src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMD.C b/src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMD.C
index e061a1a1b7..a1727fa865 100644
--- a/src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMD.C
+++ b/src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMD.C
@@ -834,8 +834,15 @@ Foam::DMDModels::STDMD::STDMD
     amps_(Zero),
     mags_(Zero),
     magsi_(Zero),
+    patches_
+    (
+        dict.getOrDefault<wordRes>
+        (
+            "patches",
+            dict.found("patch") ? wordRes(1,dict.get<word>("patch")) : wordRes()
+        )
+    ),
     fieldName_(dict.get<word>("field")),
-    patch_(dict.getOrDefault<word>("patch", word::null)),
     timeName0_(),
     minBasis_(0),
     minEval_(0),
diff --git a/src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMD.H b/src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMD.H
index f95d539b3a..006a663152 100644
--- a/src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMD.H
+++ b/src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMD.H
@@ -275,12 +275,12 @@ class STDMD
         //- Indices of 'mags'
         List<label> magsi_;
 
+        //- Names of operand patches
+        const wordRes patches_;
+
         //- Name of operand field
         const word fieldName_;
 
-        //- Name of operand patch
-        const word patch_;
-
         //- First-processed snapshot required by the mode-sorting
         //- algorithms at the final output computations (K:p. 43)
         word timeName0_;
@@ -411,7 +411,8 @@ class STDMD
                 GeoFieldType& modeRe,
                 GeoFieldType& modeIm,
                 const RMatrix& primitiveMode,
-                const label i
+                const label magi,
+                const label rowi = 0
             );
 
             //- Compute a mode for a given non-scalar-based input field
@@ -425,7 +426,8 @@ class STDMD
                 GeoFieldType& modeRe,
                 GeoFieldType& modeIm,
                 const RMatrix& primitiveMode,
-                const label i
+                const label magi,
+                const label rowi = 0
             );
 
 
diff --git a/src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMDTemplates.C b/src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMDTemplates.C
index 8dd7112150..7a68139c34 100644
--- a/src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMDTemplates.C
+++ b/src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMDTemplates.C
@@ -109,7 +109,7 @@ bool Foam::DMDModels::STDMD::calcModes()
     RxInv_.clear();
 
     label modei = 0;
-    for (const label i : magsi_)
+    for (const label magi : magsi_)
     {
         GeoFieldType modeRe
         (
@@ -147,21 +147,30 @@ bool Foam::DMDModels::STDMD::calcModes()
 
         if (modeRe.size() != 0 && !empty_)
         {
-            if (patch_.empty())
+            if (patches_.empty())
             {
                 auto& re = modeRe.primitiveFieldRef();
                 auto& im = modeIm.primitiveFieldRef();
 
-                calcMode(re, im, primitiveMode, i);
+                calcMode(re, im, primitiveMode, magi);
             }
             else
             {
-                const label patchi = mesh_.boundaryMesh().findPatchID(patch_);
+                label rowi = 0;
+                const labelList patchis
+                (
+                    mesh_.boundaryMesh().patchSet(patches_).sortedToc()
+                );
 
-                auto& re = modeRe.boundaryFieldRef()[patchi];
-                auto& im = modeIm.boundaryFieldRef()[patchi];
+                for (const label patchi : patchis)
+                {
+                    auto& re = modeRe.boundaryFieldRef()[patchi];
+                    auto& im = modeIm.boundaryFieldRef()[patchi];
 
-                calcMode(re, im, primitiveMode, i);
+                    calcMode(re, im, primitiveMode, magi, rowi);
+
+                    rowi += re.size()*pTraits<Type>::nComponents;
+                }
             }
         }
 
@@ -184,21 +193,22 @@ typename std::enable_if
     GeoFieldType& modeRe,
     GeoFieldType& modeIm,
     const RMatrix& primitiveMode,
-    const label i
+    const label magi,
+    const label rowi
 )
 {
-    const label fieldSize = modeRe.size();
+    const label szfld = modeRe.size();
 
-    for (label p = 0; p < primitiveMode.m(); ++p)
+    for (label i = rowi; i < szfld + rowi; ++i)
     {
         complex mode(Zero);
-        for (label q = 0; q < evecs_.m(); ++q)
+        for (label j = 0; j < evecs_.m(); ++j)
         {
-            mode += primitiveMode(p, q)*evecs_(q, i);
+            mode += primitiveMode(i, j)*evecs_(j, magi);
         }
-        label p1 = p%fieldSize;
-        modeRe[p1] = mode.real();
-        modeIm[p1] = mode.imag();
+        const label k = (i-rowi)%szfld;
+        modeRe[k] = mode.real();
+        modeIm[k] = mode.imag();
     }
 }
 
@@ -213,22 +223,25 @@ typename std::enable_if
     GeoFieldType& modeRe,
     GeoFieldType& modeIm,
     const RMatrix& primitiveMode,
-    const label i
+    const label magi,
+    const label rowi
 )
 {
-    const label fieldSize = modeRe.size();
+    const label szfld = modeRe.size();
+    const label szfldcmps =
+        szfld*pTraits<typename GeoFieldType::value_type>::nComponents;
 
-    for (label p = 0; p < primitiveMode.m(); ++p)
+    for (label i = rowi; i < szfldcmps + rowi; ++i)
     {
         complex mode(Zero);
-        for (label q = 0; q < evecs_.m(); ++q)
+        for (label j = 0; j < evecs_.m(); ++j)
         {
-            mode += primitiveMode(p, q)*evecs_(q, i);
+            mode += primitiveMode(i, j)*evecs_(j, magi);
         }
-        label p1 = p%fieldSize;
-        label p2 = p/fieldSize;
-        modeRe[p1][p2] = mode.real();
-        modeIm[p1][p2] = mode.imag();
+        const label k = (i-rowi)%szfld;
+        const label m = (i-rowi)/szfld;
+        modeRe[k][m] = mode.real();
+        modeIm[k][m] = mode.imag();
     }
 }
 
diff --git a/src/functionObjects/field/DMD/DMDTemplates.C b/src/functionObjects/field/DMD/DMDTemplates.C
index 3bc7016141..5b798768b1 100644
--- a/src/functionObjects/field/DMD/DMDTemplates.C
+++ b/src/functionObjects/field/DMD/DMDTemplates.C
@@ -73,35 +73,40 @@ bool Foam::functionObjects::DMD::storeSnapshot()
 
     const GeoFieldType& field = lookupObject<GeoFieldType>(fieldName_);
 
-    if (patch_.empty())
+    label rowi = nSnap_;
+    if (patches_.empty())
     {
         const label nField = field.size();
 
         for (direction dir = 0; dir < nComps; ++dir)
         {
-            z_.subColumn(0, nSnap_ + dir*nField, nField) = field.component(dir);
+            z_.subColumn(0, rowi, nField) = field.component(dir);
+            rowi += nField;
         }
     }
     else
     {
-        const label patchi = mesh_.boundaryMesh().findPatchID(patch_);
+        const labelList patchis
+        (
+            mesh_.boundaryMesh().patchSet(patches_).sortedToc()
+        );
 
-        if (patchi < 0)
+        for (const label patchi : patchis)
         {
-            FatalErrorInFunction
-                << "Cannot find patch " << patch_
-                << exit(FatalError);
-        }
+            const typename GeoFieldType::Boundary& bf = field.boundaryField();
 
-        const typename GeoFieldType::Boundary& bf = field.boundaryField();
+            const Field<typename GeoFieldType::value_type>& pbf = bf[patchi];
 
-        const Field<typename GeoFieldType::value_type>& pbf = bf[patchi];
+            const label nField = pbf.size();
 
-        const label nField = pbf.size();
-
-        for (direction dir = 0; dir < nComps; ++dir)
-        {
-            z_.subColumn(0, nSnap_ + dir*nField, nField) = pbf.component(dir);
+            if (nField > 0)
+            {
+                for (direction dir = 0; dir < nComps; ++dir)
+                {
+                    z_.subColumn(0, rowi, nField) = pbf.component(dir);
+                    rowi += nField;
+                }
+            }
         }
     }