Combined PR for number feature and io-stream wrap fix

.github/workflows/main.yml

@@ -5,7 +5,7 @@ on:
     branches: [ test ]
   pull_request:
     branches: [ test ]
-    types: [opened]
+    types: [synchronize]
 
 jobs:
   build:

.github/workflows/release_binary.yml

@@ -8,7 +8,7 @@ on:
 jobs:
   vm:
     runs-on: ubuntu-latest
-    container: quay.io/pypa/manylinux2010_x86_64    
+    container: quay.io/pypa/manylinux2014_x86_64
     steps:
       - name: build_sdist
         run: |
@@ -41,7 +41,7 @@ jobs:
         pythonpath: ["cp36-cp36m", "cp37-cp37m", "cp38-cp38", "cp39-cp39"]
 
     runs-on: ubuntu-latest
-    container: quay.io/pypa/manylinux2010_x86_64
+    container: quay.io/pypa/manylinux2014_x86_64
     #container: node:10.16-jessie
     steps:
       - name: build package
@@ -63,6 +63,7 @@ jobs:
           fi	  
           CWD=$PWD
           ls -l 
+          yum install -y zlib-devel
           python3 setup.py install
           rm -rf dist/*
           python3 setup.py bdist_wheel

.github/workflows/testrelease_binary.yml

@@ -8,7 +8,7 @@ on:
 jobs:
   vm:
     runs-on: ubuntu-latest
-    container: quay.io/pypa/manylinux2010_x86_64    
+    container: quay.io/pypa/manylinux2014_x86_64
     steps:
       - name: build_sdist
         run: |
@@ -32,7 +32,7 @@ jobs:
           python3 setup.py sdist
           ls -l dist/
           python3 -m twine upload --repository-url https://test.pypi.org/legacy/ --password ${{ secrets.TEST_PYPI_TOKEN }} --username __token__ --verbose dist/*
-          #python3 -m twine upload --password ${{ secrets.PYPI_TOKEN }} --username __token__ --verbose dist/*   	  
+          #python3 -m twine upload --password ${{ secrets.PYPI_TOKEN }} --username __token__ --verbose dist/*
   container:
     needs: vm
     strategy:
@@ -41,7 +41,7 @@ jobs:
         pythonpath: ["cp36-cp36m", "cp37-cp37m", "cp38-cp38", "cp39-cp39"]
 
     runs-on: ubuntu-latest
-    container: quay.io/pypa/manylinux2010_x86_64
+    container: quay.io/pypa/manylinux2014_x86_64
     #container: node:10.16-jessie
     steps:
       - name: build package
@@ -62,7 +62,8 @@ jobs:
              pip3 install -r requirements.txt
           fi	  
           CWD=$PWD
-          ls -l 
+          ls -l
+          yum install -y zlib-devel
           python3 setup.py install
           rm -rf dist/*
           python3 setup.py bdist_wheel

README.md

@@ -2,9 +2,9 @@
 [![badge](https://github.com/GLVis/pyglvis/blob/master/examples/ex1.svg "MFEM's Example 1")](https://mybinder.org/v2/gh/GLVis/pyglvis/HEAD?filepath=examples%2Fex1.ipynb)
 [![badge](https://github.com/GLVis/pyglvis/blob/master/examples/ex9.svg "MFEM's Example 9")](https://mybinder.org/v2/gh/GLVis/pyglvis/HEAD?filepath=examples%2Fex9.ipynb)
 
-#  PyMFEM (MFEM Python wrapper)
+#  MFEM + PyMFEM (FEM library)
 
-This package (PyMFEM) is Python wrapper for the MFEM high performance parallel finite element method library.(http://mfem.org/).
+This package provides MFEM and its Python wrapper (PyMFEM). MFEM is a high performance parallel finite element method (FEM) library (http://mfem.org/).
 
 Installer downloads a couple of external libraries and build them.
 By default, "pip install mfem" downloads and builds the serial version of MFEM and PyMFEM.

binder/requirements.txt

@@ -3,5 +3,5 @@ numpy>=1.20.0; python_version>="3.7"
 scipy>=1.5.2
 six>1.13.0
 matplotlib
-mfem
+mfem==4.2.0.15
 glvis

docs/changelog.txt

@@ -1,4 +1,16 @@
 	<<<  Change Log. >>>
+
+2021 05.19
+        * Fixed std::istream & wrapping not to return the lenght of data. This fixes
+	  the constructors such as mfem::Mesh(std::stream &)
+	* PrintGZ, SaveGZ, etc are added to the methods to save data to file. As
+	  the name suggest, they turn on compression internally.
+
+
+2021 05.11
+        * NumbaFunction, VectorNumbaFunction, and MatrixNumbaFunction is added
+	  to use Numba JITed python code for mfem function coefficients.
+
 2021 05.07
         * vtk.py and datacollection.py was loaded in mfem.ser and mfem.par namespace properly
 	* istream& wrapping was improved in the same way as ostream&

docs/manual.txt

@@ -117,6 +117,21 @@ to find more details of the MFEM library.
           Note: the argument is first tested if it can be converted using
 	        np.array(x,  dtype=float)
 
+       (EXPERIMENTAL)
+       Using Numba JIT functions for FunctionCoefficients is supported.
+       For example, to use Numba compiled function, one can use
+       GenerateCoefficient mehtod in NumberFunction class as follows.
+
+       >>> from numba import cfunc, carray
+       >>> from mfem.coefficient import scalar_sig
+       >>> sdim = mesh.SpaceDimension()
+       >>> vdim = mesh.Dimension()
+       >>> @cfunc(scalar_sig)
+           def s_func(ptx, sdim):
+           ...    return 2*ptx[0]
+       >>> c = NumbaFunction(s_func, sdim).GenerateCoefficient()
+
+
   4-3) mfem::GridFunction
        GetNodalValues(i) will perform GetNodalValue(Vector(), i) and
        return NumPy array of Vector()
@@ -295,7 +310,11 @@ to find more details of the MFEM library.
 
               # Assigning (copy) from NumPy array
 	      m = DenseMatrix(3,3)
-	      m.Assign(np.arange(9.).reshape(3,3))
+	      m.Assign(np.arange(9.).reshape(3,3))  # <-- In this case,
+	                                            #     wrapper automatically tranpose
+	                                            #     input matrix so that data is
+						    #     column order. The same is done
+						    #     in GetDataArray()
 
               x = DenseTensor(2, 3, 5)
 	      x.Assign(1) # set all 1
@@ -493,7 +512,7 @@ to find more details of the MFEM library.
             strumpack.SetOperator(Arow)
             strumpack.SetFromCommandLine()
             strumpack.Mult(B, X);
-	     
+
 -- Known issues --
 
   PyMFEM is an on-going effort. Not all MFEM functionality is properly

examples/ex4.py

@@ -6,7 +6,7 @@
 from mfem import path
 import mfem.ser as mfem
 from mfem.ser import intArray
-from os.path import expanduser, join
+from os.path import expanduser, join, dirname, exists
 import numpy as np
 from numpy import sin, array, cos
 
@@ -18,7 +18,12 @@
 kappa = np.pi*freq
 
 order = 1
+
+
 meshfile = expanduser(join(path, 'data', 'star.mesh'))
+if not exists(meshfile):
+    path = dirname(dirname(__file__))
+    meshfile = expanduser(join(path, 'data', 'star.mesh'))
 
 mesh = mfem.Mesh(meshfile, 1,1)
 dim = mesh.Dimension()

examples/numba/ex4.py

@@ -0,0 +1,122 @@
+'''
+   MFEM example 4
+
+   See c++ version in the MFEM library for more detail 
+'''
+from mfem import path
+import mfem.ser as mfem
+from mfem.ser import intArray
+from os.path import expanduser, join, dirname, exists
+import numpy as np
+from numpy import sin, array, cos
+
+set_bc = True
+static_cond = False
+hybridization = False
+
+freq = 1.0
+kappa = np.pi*freq
+
+order = 1
+
+
+meshfile = expanduser(join(path, 'data', 'star.mesh'))
+if not exists(meshfile):
+    path = dirname(dirname(__file__))
+    path = dirname(path)    
+    meshfile = expanduser(join(path, 'data', 'star.mesh'))
+
+mesh = mfem.Mesh(meshfile, 1,1)
+dim = mesh.Dimension()
+sdim = mesh.SpaceDimension()
+ndim = mesh.Dimension()
+
+from numba import cfunc, carray    
+
+@cfunc(mfem.vector_sig)
+def E_exact(p, out, sdim, vdim):
+    out_array = carray(out, (vdim, ))    
+    x = p[0]; y = p[1]
+    F0 = cos(kappa*x)*sin(kappa*y)
+    F1 = cos(kappa*y)*sin(kappa*x)
+    out_array[0] = F0
+    out_array[1] = F1
+    if vdim == 3:
+        out_array[2] = 0.0
+
+@cfunc(mfem.vector_sig)        
+def f_exact(p, out, sdim, vdim):
+    out_array = carray(out, (vdim, ))    
+    x = p[0]; y = p[1]
+    temp = 1. + 2.*kappa*kappa
+
+    F0 = temp * cos(kappa*x)*sin(kappa*y)
+    F1 = temp * cos(kappa*y)*sin(kappa*x)
+    out_array[0] = F0
+    out_array[1] = F1
+    if vdim == 3:
+        out_array[2] = 0.0
+
+ref_levels = int(np.floor(np.log(25000./mesh.GetNE())/np.log(2.)/dim))
+for x in range(ref_levels):
+   mesh.UniformRefinement();
+
+fec = mfem.RT_FECollection(order-1, dim)
+fespace = mfem.FiniteElementSpace(mesh, fec)
+
+print("Number of finite element unknows : " + str(fespace.GetTrueVSize()))
+
+ess_tdof_list = intArray();
+if mesh.bdr_attributes.Size():
+    ess_bdr = intArray(mesh.bdr_attributes.Max())
+    if set_bc: ess_bdr.Assign(1)
+    else: ess_bdr.Assign(0)
+    fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+
+b = mfem.LinearForm(fespace);
+f = mfem.VectorNumbaFunction(f_exact, sdim, ndim).GenerateCoefficient()
+dd = mfem.VectorFEDomainLFIntegrator(f)
+b.AddDomainIntegrator(dd)
+b.Assemble()
+
+x = mfem.GridFunction(fespace)
+F = mfem.VectorNumbaFunction(E_exact, sdim, ndim).GenerateCoefficient()
+x.ProjectCoefficient(F);
+
+alpha = mfem.ConstantCoefficient(1.0);
+beta = mfem.ConstantCoefficient(1.0);
+a = mfem.BilinearForm(fespace);
+a.AddDomainIntegrator(mfem.DivDivIntegrator(alpha));
+a.AddDomainIntegrator(mfem.VectorFEMassIntegrator(beta));
+
+
+if (static_cond):  a.EnableStaticCondensation()
+elif (hybridization):
+    hfec = mfem.DG_Interface_FECollection(order-1, dim)
+    hfes = mfem.FiniteElementSpace(mesh, hfec)
+    a.EnableHybridization(hfes, mfem.NormalTraceJumpIntegrator(),
+                          ess_tdof_list)
+a.Assemble();
+
+
+A = mfem.OperatorPtr()
+B = mfem.Vector()
+X = mfem.Vector()
+a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+## Here, original version calls hegith, which is not
+## defined in the header...!?
+print("Size of linear system: " + str(A.Height())) 
+
+# 10. Solve
+AA = mfem.OperatorHandle2SparseMatrix(A)
+M = mfem.GSSmoother(AA)
+X.Print("x")
+mfem.PCG(AA, M, B, X, 1, 10000, 1e-20, 0.0);
+
+# 11. Recover the solution as a finite element grid function.
+a.RecoverFEMSolution(X, b, x);
+
+print("|| F_h - F ||_{L^2} = " + str(x.ComputeL2Error(F)))
+
+mesh.Print('refined.mesh', 8)
+x.Save('sol.gf', 8)

mfem/__init__.py

@@ -20,6 +20,4 @@ def debug_print(message):
 
     print(message)
 
-__version__ = '4.2.0.14'
-
-
+__version__ = '4.2.0.16'

mfem/_par/array.py

@@ -348,6 +348,16 @@ def Print(self, *args):
         return _array.intArray_Print(self, *args)
     Print = _swig_new_instance_method(_array.intArray_Print)
 
+    def PrintGZ(self, file, precision=8):
+        r"""PrintGZ(intArray self, char const * file, int precision=8)"""
+        return _array.intArray_PrintGZ(self, file, precision)
+    PrintGZ = _swig_new_instance_method(_array.intArray_PrintGZ)
+
+    def SaveGZ(self, file, precision=8):
+        r"""SaveGZ(intArray self, char const * file, int precision=8)"""
+        return _array.intArray_SaveGZ(self, file, precision)
+    SaveGZ = _swig_new_instance_method(_array.intArray_SaveGZ)
+
     def Save(self, *args):
         r"""
         Save(intArray self, std::ostream & out, int fmt=0)
@@ -643,6 +653,16 @@ def Print(self, *args):
         return _array.doubleArray_Print(self, *args)
     Print = _swig_new_instance_method(_array.doubleArray_Print)
 
+    def PrintGZ(self, file, precision=8):
+        r"""PrintGZ(doubleArray self, char const * file, int precision=8)"""
+        return _array.doubleArray_PrintGZ(self, file, precision)
+    PrintGZ = _swig_new_instance_method(_array.doubleArray_PrintGZ)
+
+    def SaveGZ(self, file, precision=8):
+        r"""SaveGZ(doubleArray self, char const * file, int precision=8)"""
+        return _array.doubleArray_SaveGZ(self, file, precision)
+    SaveGZ = _swig_new_instance_method(_array.doubleArray_SaveGZ)
+
     def Save(self, *args):
         r"""
         Save(doubleArray self, std::ostream & out, int fmt=0)