Merge pull request #101 from mfem/test

Fix memory leak issue in DenseMatrix::Assign

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==4.3.0.0
+mfem==4.3.0.1
 glvis

docs/changelog.txt

@@ -1,13 +1,18 @@
 	<<<  Change Log. >>>
+2021 08.27
+        * memory leak
+        * DenseTensor::Assign now accept numpy 3D array
+
 2021 08.08
         * added ex0p, ex20p, ex21p, ex22p  ex24p, ex26p, ex27p, ex28p
-	* memory handling of complex_fem for parallel version is fixed.
+        * memory handling of complex_fem for parallel version is fixed.
 
 2021 08.06
         * ex27 is added
 	* ex26 is added
 	* array((int *, size)) and array((double *, size)) is wrapped. This allows
 	  for changing the contents of exisint Array<T>
+
 2021 08.04
 	* added ex24. This one shows partial assembly and numba coefficient
 

docs/manual.txt

@@ -43,6 +43,11 @@ to find more details of the MFEM library.
 	>>> print(x[0])   # access to element
 	>>> x[3] = 5      # set element
 
+     mfem.intArray provides an access to int *. For this, we need to
+     pass (int *, length) as a list/tuple.
+        >>> parts = mesh.GeneratePartitioning(3)
+        >>> mfem.intArray([parts, mesh.GetNE()])   <-- int * and length 
+
 
   4-1) mfem::Vector
      Constructor using python list or NumPy array:
@@ -365,8 +370,23 @@ to find more details of the MFEM library.
 	      x[1,2,3]=3
 
               # NOTE the index order is different in the NumPy array returned
-	      # by DenseTensor::GetDataArray
+	      # by DenseTensor::GetDataArray. This is because, we want 
+	      # the major column access work in the same way in C and Python
+	      # x[3].GetDataArray() == x.GetDataArray()[3]
+
               x[1 2,3] = x.GetDataArray()[3,1,2]
+
+              # tensor assignment behaves as expected. 
+              t = np.arange(30).reshape(2,3,5)
+	      x.Assign(t) # tensor assignment
+
+              # because of the same reason as DenseMatrix, this operation changes
+	      # internall array element order
+	      t.flatten() ---> [0, 1, 2, 3, 4...]
+	      print(mfem.Vector(m.Data(), 30).GetDataArray()) ---> [0, 15, 5,,,,]
+
+
+
 
 	 Note that internally MFEM DenseMatrix is column major, while NumPy
 	 is row major. Therefore,

examples/ex18.py

@@ -201,7 +201,7 @@
     sol_name = "vortex-" + str(k) + "-final.gf"
     uk.Save(sol_name, 8)
 
-print("done")
+print(" done")
 # 10. Compute the L2 solution error summed for all components.
 if (t_final == 2.0):
     error = sol.ComputeLpError(2., u0)

examples/ex18_common.py

@@ -56,13 +56,20 @@ def GetFlux(self, x, flux):
         state = self.state
         dof = self.flux.SizeI()
         dim = self.flux.SizeJ()
+
+        flux_data = []
         for i in range(dof):
             for k in range(num_equation):
                 self.state[k] = x[i, k]
             ComputeFlux(state, dim, self.f)
-            for d in range(dim):
-                for k in range(num_equation):
-                    flux[i, d, k] = self.f[k, d]
+
+            flux_data.append(self.f.GetDataArray().transpose().copy())
+            #flux[i].Print()
+            #print(self.f.GetDataArray())
+            #for d in range(dim):
+            #    for k in range(num_equation):
+            #        flux[i, d, k] = self.f[k, d]
+        flux.Assign(np.stack(flux_data))
         mcs = ComputeMaxCharSpeed(state, dim)
         if (mcs > max_char_speed):
             globals()['max_char_speed'] = mcs
@@ -282,7 +289,7 @@ def StateIsPhysical(state, dim):
     specific_heat_ratio = globals()["specific_heat_ratio"]
 
     den = state[0]
-    den_vel = state[1:1+dim].GetDataArray()
+    den_vel = state.GetDataArray()[1:1+dim]
     den_energy = state[1 + dim]
 
     if (den < 0):
@@ -361,7 +368,7 @@ def EvalValue(self, x):
 
 def ComputePressure(state, dim):
     den = state[0]
-    den_vel = state[1:1+dim].GetDataArray()
+    den_vel = state.GetDataArray()[1:1+dim]
     den_energy = state[1 + dim]
 
     specific_heat_ratio = globals()["specific_heat_ratio"]
@@ -370,14 +377,13 @@ def ComputePressure(state, dim):
 
     return pres
 
-
 def ComputeFlux(state, dim, flux):
     from mpi4py import MPI
     num_procs = MPI.COMM_WORLD.size
     myid = MPI.COMM_WORLD.rank
 
     den = state[0]
-    den_vel = state[1:1+dim].GetDataArray()
+    den_vel = state.GetDataArray()[1:1+dim]
     den_energy = state[1 + dim]
 
     assert StateIsPhysical(state, dim), ""
@@ -391,14 +397,6 @@ def ComputeFlux(state, dim, flux):
     for d in range(dim):
         fluxA[1+d, d] += pres
 
-    '''
-    for d in range(dim):
-       flux[0, d] = den_vel[d]
-       for i in range(dim):       
-          flux[1+i, d] = den_vel[i] * den_vel[d] / den;
-       flux[1+d, d] += pres;
-    '''
-
     H = (den_energy + pres) / den
     flux.GetDataArray()[1+dim, :] = den_vel * H
 
@@ -410,7 +408,7 @@ def ComputeFluxDotN(state, nor, fluxN):
     fluxN = fluxN.GetDataArray()
 
     den = state[0]
-    den_vel = state[1:1+dim].GetDataArray()
+    den_vel = state.GetDataArray()[1:1+dim]    
     den_energy = state[1 + dim]
 
     assert StateIsPhysical(state, dim), ""
@@ -430,7 +428,7 @@ def ComputeMaxCharSpeed(state, dim):
     specific_heat_ratio = globals()["specific_heat_ratio"]
 
     den = state[0]
-    den_vel = state[1:1+dim].GetDataArray()
+    den_vel = state.GetDataArray()[1:1+dim]    
     den_energy = state[1 + dim]
 
     den_vel2 = np.sum(den_vel**2)

mfem/__init__.py

@@ -20,4 +20,4 @@ def debug_print(message):
 
     print(message)
 
-__version__ = '4.3.0.0'
+__version__ = '4.3.0.1'

mfem/_par/densemat.i

@@ -133,7 +133,7 @@ def __getitem__(self, *args):
     }
     int ndim = PyArray_NDIM(numpymat);
     if (ndim != 2){
-      PyErr_SetString(PyExc_ValueError, "Input data NDIM must be two");
+      PyErr_SetString(PyExc_ValueError, "Input data NDIM must be 2");
       return ;
     }
     npy_intp *shape = PyArray_DIMS(numpymat);
@@ -142,10 +142,13 @@ def __getitem__(self, *args):
       PyErr_SetString(PyExc_ValueError, "input data length does not match");
       return ;
     }
-    PyArrayObject * tmp = 
-      PyArray_GETCONTIGUOUS((PyArrayObject *)PyArray_Transpose((PyArrayObject *)numpymat, NULL));
-    (* self) = (double *) PyArray_DATA(tmp);
-    Py_XDECREF(tmp);
+    PyObject * tmp1 = 
+       PyArray_Transpose((PyArrayObject *)numpymat, NULL);
+    PyArrayObject * tmp2 = 
+      PyArray_GETCONTIGUOUS((PyArrayObject *)tmp1);
+    (* self) = (double *) PyArray_DATA(tmp2);
+    Py_XDECREF(tmp1);
+    Py_XDECREF(tmp2);    
   }
 
   const double __getitem__(const int i, const int j) const{
@@ -157,14 +160,53 @@ def __getitem__(self, *args):
   PyObject* GetDataArray(void) const{
      double * A = self->Data();    
      npy_intp dims[] = {self->Width(), self->Height()};
-     return  PyArray_Transpose((PyArrayObject *)PyArray_SimpleNewFromData(2, dims, NPY_DOUBLE, A), NULL);
+     PyObject *tmp1 = PyArray_SimpleNewFromData(2, dims, NPY_DOUBLE, A);
+     PyObject *ret = PyArray_Transpose((PyArrayObject *)tmp1, NULL);
+     Py_XDECREF(tmp1);
+     return ret;
+       //return  PyArray_Transpose((PyArrayObject *)PyArray_SimpleNewFromData(2, dims, NPY_DOUBLE, A), NULL);
   }
 };
 
 %extend mfem::DenseTensor {
   void Assign(const double c) {
     (* self) = c;
   }
+  void Assign(const mfem::DenseTensor &m) {
+    (* self) = m;
+  }
+  void Assign(PyObject* numpymat) {
+    /* note that these error does not raise error in python
+       type check is actually done in wrapper layer */
+    if (!PyArray_Check(numpymat)){
+       PyErr_SetString(PyExc_ValueError, "Input data must be ndarray");
+       return;
+    }
+    int typ = PyArray_TYPE(numpymat);
+    if (typ != NPY_DOUBLE){
+        PyErr_SetString(PyExc_ValueError, "Input data must be float64");
+	return;
+    }
+    int ndim = PyArray_NDIM(numpymat);
+    if (ndim != 3){
+      PyErr_SetString(PyExc_ValueError, "Input data NDIM must be 3");
+      return ;
+    }
+    npy_intp *shape = PyArray_DIMS(numpymat);    
+    int len = self->SizeI()*self->SizeJ()*self->SizeK();
+    if (shape[2]*shape[1]*shape[0] != len){    
+      PyErr_SetString(PyExc_ValueError, "input data length does not match");
+      return ;
+    }
+
+    for (int i=0; i < self->SizeI(); i++){
+       for (int j=0; j < self->SizeJ(); j++){
+          for (int k=0; k < self->SizeK(); k++){      
+	    (* self)(i, j, k) = *(double *) PyArray_GETPTR3((PyArrayObject *)numpymat, i, j, k);
+	}
+      }
+    }
+  }
   const double __getitem__(const int i, const int j, const int k) const{
     return (* self)(i, j, k);
   }
@@ -180,8 +222,9 @@ def __getitem__(self, *args):
      npy_intp dims[] = {self->SizeK(), self->SizeJ(), self->SizeI()};
      PyObject * obj = PyArray_SimpleNewFromData(3, dims, NPY_DOUBLE, A);
      //obj = PyArray_SwapAxes((PyArrayObject *)obj, 0, 2);
-     obj = PyArray_SwapAxes((PyArrayObject *)obj, 1, 2);
-     return obj;
+     PyObject * ret = PyArray_SwapAxes((PyArrayObject *)obj, 1, 2);
+     Py_XDECREF(obj);     
+     return ret;
   }
 };
 

mfem/_par/densemat.py

@@ -1134,9 +1134,13 @@ def Swap(self, t):
     Swap = _swig_new_instance_method(_densemat.DenseTensor_Swap)
     __swig_destroy__ = _densemat.delete_DenseTensor
 
-    def Assign(self, c):
-        r"""Assign(DenseTensor self, double const c)"""
-        val = _densemat.DenseTensor_Assign(self, c)
+    def Assign(self, *args):
+        r"""
+        Assign(DenseTensor self, double const c)
+        Assign(DenseTensor self, DenseTensor m)
+        Assign(DenseTensor self, PyObject * numpymat)
+        """
+        val = _densemat.DenseTensor_Assign(self, *args)
 
         return self