Added an API to read OF fields.

Author: friedenhe

pyofm/__init__.py

@@ -1,3 +1,3 @@
-__version__ = '1.2.0'
+__version__ = '1.2.1'
 
 from .pyOFM import PYOFM

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)

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);
+        }
+    }
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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)
+        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)

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",