From 24f511be90f1510aca834adf438897dd10fdee85 Mon Sep 17 00:00:00 2001
From: Ping He <phe@iastate.edu>
Date: Sun, 16 Jun 2024 16:57:32 -0500
Subject: [PATCH] Added an API to read OF fields.

---
 pyofm/__init__.py |  2 +-
 pyofm/pyOFM.py    | 14 ++++++++++++
 src/OFMesh.H      | 56 +++++++++++++++++++++++++++++++++++++++++++++--
 src/pyOFMesh.pyx  | 17 ++++++++++++--
 src/setup.py      |  3 +++
 5 files changed, 87 insertions(+), 5 deletions(-)

diff --git a/pyofm/__init__.py b/pyofm/__init__.py
index 656569d..f714086 100755
--- a/pyofm/__init__.py
+++ b/pyofm/__init__.py
@@ -1,3 +1,3 @@
-__version__ = '1.2.0'
+__version__ = '1.2.1'
 
 from .pyOFM import PYOFM
diff --git a/pyofm/pyOFM.py b/pyofm/pyOFM.py
index 8ac02c9..d0951a3 100755
--- a/pyofm/pyOFM.py
+++ b/pyofm/pyOFM.py
@@ -306,3 +306,17 @@ def writeCompression(self):
                 if "on" in cols[1]:
                     writecompress = True
         return writecompress
+
+    def readField(self, fieldName, fieldType, timeName, field):
+        """
+        Read OpenFoam field and return the internal field as an array
+
+        Inputs:
+            fieldName: name of the field to read
+            fieldType: can be either volScalarField or volVectorField
+            timeName: the time folder name to read, e.g., "0" or "1000"
+        Output:
+            field: an np array to save the field
+        """
+
+        self.ofMesh.readField(fieldName, fieldType, timeName, field)
diff --git a/src/OFMesh.H b/src/OFMesh.H
index 5ab2ebc..583ac2e 100755
--- a/src/OFMesh.H
+++ b/src/OFMesh.H
@@ -35,12 +35,12 @@ private:
     /// all the arguments
     char* argsAll_;
 
-    autoPtr<argList> argsPtr_;
-
     autoPtr<fvMesh> meshPtr_;
 
     autoPtr<Time> runTimePtr_;
 
+    autoPtr<argList> argsPtr_;
+
     label nLocalPoints_;
 
     label nLocalCells_;
@@ -171,6 +171,58 @@ public:
         const UList<label>& pFaceCells = meshPtr_().boundaryMesh()[patchI].faceCells();
         return pFaceCells[faceI];
     }
+
+    void readField(
+        const word fieldName,
+        const word fieldType,
+        const word timeName,
+        double* field)
+    {
+        if (fieldType == "volScalarField")
+        {
+            volScalarField state(
+                IOobject(
+                    fieldName,
+                    timeName,
+                    meshPtr_(),
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE,
+                    false),
+                meshPtr_());
+
+            forAll(state, cellI)
+            {
+                field[cellI] = state[cellI];
+            }
+        }
+        else if (fieldType == "volVectorField")
+        {
+            volVectorField state(
+                IOobject(
+                    fieldName,
+                    timeName,
+                    meshPtr_(),
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE,
+                    false),
+                meshPtr_());
+
+            label counterI = 0;
+            forAll(state, cellI)
+            {
+                for (label i = 0; i < 3; i++)
+                {
+                    field[counterI] = state[cellI][i];
+                    counterI++;
+                }
+            }
+        }
+        else
+        {
+            FatalErrorIn("readField") << "fieldType " << fieldType
+                                      << " not supported! Options are volScalariField or volVectorField" << abort(FatalError);
+        }
+    }
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/pyOFMesh.pyx b/src/pyOFMesh.pyx
index 8f7911a..98e3b46 100755
--- a/src/pyOFMesh.pyx
+++ b/src/pyOFMesh.pyx
@@ -14,7 +14,7 @@
 
 '''
 
-import numpy as np
+cimport numpy as np
 
 from libcpp.string cimport string
 
@@ -41,6 +41,7 @@ cdef extern from "OFMesh.H" namespace "Foam":
         int getLocalFaceOwner(int)
         int getLocalBoundaryFaceOwner(int,int)
         int getLocalFaceNeighbour(int)
+        void readField(char* , char *, char *, double *)
 
 # create python wrappers that call cpp functions
 cdef class pyOFMesh:
@@ -132,4 +133,16 @@ cdef class pyOFMesh:
         return self._thisptr.getLocalBoundaryFaceOwner(patchI,faceI)
     
     def getLocalFaceNeighbour(self, faceI):
-        return self._thisptr.getLocalFaceNeighbour(faceI)
\ No newline at end of file
+        return self._thisptr.getLocalFaceNeighbour(faceI)
+    
+    def readField(self, fieldName, fieldType, timeName, np.ndarray[double, ndim=1, mode="c"] field):
+        if fieldType == "volScalarField":
+            assert len(field) == self.getNLocalCells(), "invalid array size!"
+        elif fieldType == "volVectorField":
+            assert len(field) == self.getNLocalCells() * 3, "invalid array size!"
+        else:
+            print("fieldType invalid!")
+            exit(1)
+        
+        cdef double *field_data = <double*>field.data
+        self._thisptr.readField(fieldName.encode(), fieldType.encode(), timeName.encode(), field_data)
diff --git a/src/setup.py b/src/setup.py
index c88cad7..2074f61 100755
--- a/src/setup.py
+++ b/src/setup.py
@@ -15,6 +15,7 @@
 from distutils.extension import Extension
 from Cython.Build import cythonize
 import os
+import numpy
 
 libName = "pyOFMesh"
 
@@ -47,6 +48,7 @@
             os.getenv("FOAM_SRC") + "/codipack/include",
             os.getenv("FOAM_SRC") + "/medipack/include",
             os.getenv("FOAM_SRC") + "/medipack/src",
+            numpy.get_include(),
         ],
         # These are from Make/options:EXE_LIBS
         libraries=["meshTools", "finiteVolume"],
@@ -56,6 +58,7 @@
         extra_compile_args=[
             # "-DFULLDEBUG -g -O0", # this is for debugging
             "-std=c++11",
+            "-Wno-deprecated-copy",
             "-m64",
             "-DOPENFOAM_PLUS=1812",
             "-Dlinux64",