diff --git a/Numeric/fullMatrix.cpp b/Numeric/fullMatrix.cpp
index e95c99eb27259cb70be9ee934ad69cd6e2e21fe6..c5197d2dd62b0daae67811b9d8e42712fc427b2c 100644
--- a/Numeric/fullMatrix.cpp
+++ b/Numeric/fullMatrix.cpp
@@ -66,6 +66,16 @@ void fullVector<std::complex<double> >::setAll(const fullVector<std::complex<dou
   F77NAME(zcopy)(&_r, m._data, &stride, _data, &stride);
 }
 
+template<>
+void fullMatrix<int>::setAll(const fullMatrix<int> &m)
+{
+  if (_r != m._r || _c != m._c )
+    Msg::Fatal("fullMatrix size does not match");
+  int N = _r * _c;
+  for (int i = 0; i < N; ++i)
+    _data[i] = m._data[i];
+}
+
 template<>
 void fullMatrix<double>::setAll(const fullMatrix<double> &m)
 {
@@ -213,6 +223,13 @@ void fullMatrix<double>::multWithATranspose(const fullVector<double> &x, double
                  &beta, y._data, &INCY);
 
 }
+
+
+template<>
+void fullMatrix<int>::scale(const double s)
+{
+  for(int i = 0; i < _r * _c; ++i) _data[i] *= s;
+}
 #endif
 
 
@@ -585,3 +602,41 @@ bool fullMatrix<double>::invert(fullMatrix<double> &result) const
 }
 
 #endif
+
+template<>
+void fullMatrix<double>::print(const std::string name, const std::string format) const
+{
+  std::string rformat = (format == "") ? "%12.5E " : format;
+  int ni = size1();
+  int nj = size2();
+  printf("double %s [ %d ][ %d ]= { \n", name.c_str(),ni,nj);
+  for(int I = 0; I < ni; I++){
+    printf("{  ");
+    for(int J = 0; J < nj; J++){
+      printf(rformat.c_str(), (*this)(I, J));
+      if (J!=nj-1)printf(",");
+    }
+    if (I!=ni-1)printf("},\n");
+    else printf("}\n");
+  }
+  printf("};\n");
+}
+
+template<>
+void fullMatrix<int>::print(const std::string name, const std::string format) const
+{
+  std::string rformat = (format == "") ? "%12d " : format;
+  int ni = size1();
+  int nj = size2();
+  printf("int %s [ %d ][ %d ]= { \n", name.c_str(),ni,nj);
+  for(int I = 0; I < ni; I++){
+    printf("{  ");
+    for(int J = 0; J < nj; J++){
+      printf(rformat.c_str(), (*this)(I, J));
+      if (J!=nj-1)printf(",");
+    }
+    if (I!=ni-1)printf("},\n");
+    else printf("}\n");
+  }
+  printf("};\n");
+}
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index 2bfd6f767caf849a63aa99e29a640dba919a0f19..f566a1fb6639a931fa56ca9c53b10b8708ab9d3c 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -430,7 +430,7 @@ class fullMatrix
     if (init0)
       setAll(scalar(0.));
   }
-  fullMatrix(int r, int c, double *data)
+  fullMatrix(int r, int c, scalar *data)
     : _r(r), _c(c), _data(data), _own_data(false)
   {
     setAll(scalar(0.));
@@ -526,7 +526,7 @@ class fullMatrix
     _own_data = false;
     _data = original._data + c_start * _r;
   }
-  void setAsProxy(double *data, int r, int c)
+  void setAsProxy(scalar *data, int r, int c)
   {
     if(_data && _own_data)
       delete [] _data;
@@ -774,22 +774,7 @@ class fullMatrix
 #endif
   ;
 
-  void print(const std::string name = "", const std::string format = "%12.5E ") const
-  {
-    int ni = size1();
-    int nj = size2();
-    printf("double %s [ %d ][ %d ]= { \n", name.c_str(),ni,nj);
-    for(int I = 0; I < ni; I++){
-      printf("{  ");
-      for(int J = 0; J < nj; J++){
-        printf(format.c_str(), (*this)(I, J));
-	if (J!=nj-1)printf(",");
-      }
-      if (I!=ni-1)printf("},\n");
-      else printf("}\n");
-    }
-    printf("};\n");
-  }
+  void print(const std::string name = "", const std::string format = "") const;
 
   void binarySave (FILE *f) const
   {
diff --git a/wrappers/gmshpy/gmshNumeric.i b/wrappers/gmshpy/gmshNumeric.i
index e58e8f7795193f8689736a9cf67fbb9278eda1ff..3ef62e3bedb2f9e8cf56c3815dc1e8c0c65fd41f 100644
--- a/wrappers/gmshpy/gmshNumeric.i
+++ b/wrappers/gmshpy/gmshNumeric.i
@@ -28,10 +28,26 @@
 %include "JacobianBasis.h"
 %ignore fullMatrix<double>::operator()(int, int);
 %ignore fullVector<double>::operator()(int);
+%ignore fullMatrix<int>::operator()(int, int);
+%ignore fullMatrix<int>::eig;
+%ignore fullMatrix<int>::luSubstitute;
+%ignore fullMatrix<int>::invert;
+%ignore fullMatrix<int>::invertInPlace;
+%ignore fullMatrix<int>::axpy;
+%ignore fullMatrix<int>::gemm;
+%ignore fullMatrix<int>::mult;
+%ignore fullMatrix<int>::luSolve;
+%ignore fullMatrix<int>::luFactor;
+%ignore fullMatrix<int>::svd;
+%ignore fullMatrix<int>::multAddy;
+%ignore fullMatrix<int>::multWithATranspose;
+%ignore fullMatrix<int>::multOnBlock;
+%ignore fullMatrix<int>::determinant;
 %include "fullMatrix.h"
 %include "simpleFunction.h"
 %template(fullMatrixDouble) fullMatrix<double>;
 %template(fullVectorDouble) fullVector<double>;
+%template(fullMatrixInt) fullMatrix<int>;
 %template(simpleFunctionDouble) simpleFunction<double>;
 %include "simpleFunctionPython.h"
 %include "nodalBasis.h"
diff --git a/wrappers/gmshpy/gmshtypemaps.i b/wrappers/gmshpy/gmshtypemaps.i
index 76acb7ccab9247089feced450b6dab8a896e6609..e903f71b3b7f06dfa9db7a57945cd1aec63183c1 100644
--- a/wrappers/gmshpy/gmshtypemaps.i
+++ b/wrappers/gmshpy/gmshtypemaps.i
@@ -61,7 +61,7 @@
       return fm;
     %#ifdef HAVE_NUMPY
     PyArrayObject *array = (PyArrayObject*)obj;
-    if (PyArray_Check(array) && PyArray_ISFARRAY_RO(array) && PyArray_NDIM(array) == 2) {
+    if (PyArray_Check(array) && PyArray_ISFARRAY_RO(array) && PyArray_NDIM(array) == 2 && PyArray_DESCR(array)->kind == 'f') {
       newMatrix = true;
       return new fullMatrix<double>((double*)PyArray_DATA(array), PyArray_DIM(array, 0), PyArray_DIM(array, 1));
     }
@@ -75,7 +75,7 @@
       newMatrix = true;
     return fm;
   }
-
+  
   fullMatrix<double> *objToFullMatrixRW(PyObject *obj, bool &newMatrix)
   {
     fullMatrix<double> *fm = NULL;
@@ -84,7 +84,7 @@
       return fm;
     %#ifdef HAVE_NUMPY
     PyArrayObject *array = (PyArrayObject*)obj;
-    if (PyArray_Check(array) && PyArray_ISFARRAY(array) && PyArray_NDIM(array) == 2) {
+    if (PyArray_Check(array) && PyArray_ISFARRAY(array) && PyArray_NDIM(array) == 2 && PyArray_DESCR(array)->kind == 'f') {
       newMatrix = true;
       return new fullMatrix<double>((double*)PyArray_DATA(array), PyArray_DIM(array, 0), PyArray_DIM(array, 1));
     }
@@ -152,6 +152,154 @@
   }
   %#endif
 }
+%fragment("fullMatrixConversionInt", "header", fragment="fullMatrixConversionInit") {
+  %#undef HAVE_DLOPEN
+  %#include "fullMatrix.h"
+  %#include "GmshConfig.h"
+  %#ifdef HAVE_NUMPY
+  %#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
+  %#include <numpy/arrayobject.h>
+  %#if NPY_API_VERSION < 0x00000007
+  %#define NPY_ARRAY_FARRAY NPY_FARRAY
+  %#define NPY_ARRAY_ALIGNED NPY_ALIGNED
+  %#define NPY_ARRAY_F_CONTIGUOUS NPY_F_CONTIGUOUS
+  %#define NPY_ARRAY_WRITEABLE NPY_WRITEABLE
+  %#define PyArray_SetBaseObject(a, b) PyArray_BASE(a) = b
+  %#endif
+  %#endif
+  fullMatrix<int> *pySequenceToFullMatrixInt(PyObject *o) {
+    fullMatrix<int> *fm = NULL;
+    if (!PySequence_Check(o))
+      return NULL;
+    long int nRow = PySequence_Length(o);
+    for (int i = 0; i < PySequence_Length(o); ++i) {
+      PyObject *l = PySequence_GetItem(o, i);
+      if (!PySequence_Check(l)){
+        if (fm != NULL)
+          delete fm;
+        return NULL;
+      }
+      long int nCol = PySequence_Length(l);
+      if (i == 0)
+        fm = new fullMatrix<int>(nRow, nCol);
+      else if (nCol != fm->size2()) {
+        delete fm;
+        return NULL;
+      }
+      for (int j = 0; j < nCol; ++j) {
+        PyObject *v = PySequence_GetItem(l, j);
+        if (PyLong_Check(v)) {
+          (*fm)(i, j) = (int) PyLong_AsLong(v);
+        }
+        else {
+          delete fm;
+          return NULL;
+        }
+      }
+    }
+    return fm;
+  }
+
+  fullMatrix<int> *objToFullMatrixIntRO(PyObject *obj, bool &newMatrix, PyObject *&tmpObject)
+  {
+    fullMatrix<int> *fm = NULL;
+    SWIG_ConvertPtr(obj,(void **) &fm, SWIGTYPE_p_fullMatrixT_int_t, 1);
+    if (fm)
+      return fm;
+    %#ifdef HAVE_NUMPY
+    PyArrayObject *array = (PyArrayObject*)obj;
+    if (PyArray_Check(array) && PyArray_ISFARRAY_RO(array) && PyArray_NDIM(array) == 2 && PyArray_DESCR(array)->kind == 'i') {
+      newMatrix = true;
+      return new fullMatrix<int>((int*)PyArray_DATA(array), PyArray_DIM(array, 0), PyArray_DIM(array, 1));
+    }
+    if ((tmpObject = PyArray_FROMANY(obj, NPY_INT, 2, 2, NPY_ARRAY_FARRAY))) {
+      array = (PyArrayObject*)tmpObject;
+      newMatrix = true;
+      return new fullMatrix<int>((int*)PyArray_DATA(array), PyArray_DIM(array, 0), PyArray_DIM(array, 1));
+    }
+    %#endif
+    if ((fm = pySequenceToFullMatrixInt(obj)))
+      newMatrix = true;
+    return fm;
+  }
+
+  fullMatrix<int> *objToFullMatrixIntRW(PyObject *obj, bool &newMatrix)
+  {
+    fullMatrix<int> *fm = NULL;
+    SWIG_ConvertPtr(obj,(void **) &fm, SWIGTYPE_p_fullMatrixT_int_t, 1);
+    if (fm)
+      return fm;
+    %#ifdef HAVE_NUMPY
+    PyArrayObject *array = (PyArrayObject*)obj;
+    if (PyArray_Check(array) && PyArray_ISFARRAY(array) && PyArray_NDIM(array) == 2 && PyArray_DESCR(array)->kind == 'i') {
+      newMatrix = true;
+      return new fullMatrix<int>((int*)PyArray_DATA(array), PyArray_DIM(array, 0), PyArray_DIM(array, 1));
+    }
+    %#endif
+    return NULL;
+  }
+
+  %#ifdef HAVE_NUMPY
+  %#if PY_MAJOR_VERSION==2 && PY_MINOR_VERSION==6
+  static void deleteCapsuleArrayInt(void *capsule)
+  {
+    delete [](double*)(capsule);
+  }
+  %#else
+  static void deleteCapsuleArrayInt(PyObject *capsule)
+  {
+    delete [](double*)PyCapsule_GetPointer(capsule, NULL);
+  }
+  %#endif
+
+  PyObject *fullMatrix2PyArray(fullMatrix<int> &fm)
+  {
+    npy_intp dims[2] = {fm.size1(), fm.size2()};
+    int *data = &fm.operator()(0, 0);
+    /*PyObject *array = PyArray_New(&PyArray_Type, 2, dims, NPY_DOUBLE, NULL, NULL, 0, NPY_ARRAY_F_CONTIGUOUS, NULL);
+    // copy data
+    memcpy((void*)PyArray_DATA(array), data, dims[0] * dims[1] * sizeof(int));*/
+    // do not copy data
+    PyObject *array = PyArray_New(&PyArray_Type, 2, dims, NPY_INT, NULL, (void*)data, 0, NPY_ARRAY_F_CONTIGUOUS, NULL);
+    PyArray_UpdateFlags((PyArrayObject*)array, NPY_ARRAY_ALIGNED);
+    if (fm.getOwnData()) {
+      fm.setOwnData(false);
+      %#if PY_MAJOR_VERSION==2 && PY_MINOR_VERSION==6
+      PyObject *capsule = PyCObject_FromVoidPtr((void*)data, deleteCapsuleArrayInt);
+      %#else
+      PyObject *capsule = PyCapsule_New((void*) data, NULL, deleteCapsuleArrayInt);
+      %#endif
+      PyArray_SetBaseObject((PyArrayObject*)array, capsule);
+    }
+    return array;
+  }
+  PyObject *fullMatrix2PyArrayConst(const fullMatrix<int> &fm)
+  {
+    npy_intp dims[2] = {fm.size1(), fm.size2()};
+    int *data = (int*)fm.getDataPtr();
+    // do not copy data
+    PyObject *array = PyArray_New(&PyArray_Type, 2, dims, NPY_INT, NULL, (void*)data, 0, NPY_ARRAY_F_CONTIGUOUS, NULL);
+    PyArray_UpdateFlags((PyArrayObject*)array, NPY_ARRAY_ALIGNED);
+    return array;
+  }
+  PyObject *fullMatrix2PyArrayProxy(const fullMatrix<int> &fm)
+  {
+    npy_intp dims[2] = {fm.size1(), fm.size2()};
+    int *data = (int*)fm.getDataPtr();
+    PyObject *array = PyArray_New(&PyArray_Type, 2, dims, NPY_INT, NULL, (void*)data, 0, NPY_ARRAY_F_CONTIGUOUS, NULL);
+    PyArray_UpdateFlags((PyArrayObject*)array, NPY_ARRAY_ALIGNED |NPY_ARRAY_WRITEABLE);
+    return array;
+  }
+  PyObject *fullVector2PyArrayProxy(fullVector<int> &fv)
+  {
+    npy_intp dims[1] = {fv.size()};
+    int *data = &fv.operator()(0);
+    PyObject *array = PyArray_New(&PyArray_Type, 1, dims, NPY_INT, NULL, (void*)data, 0, NPY_ARRAY_F_CONTIGUOUS, NULL);
+    PyArray_UpdateFlags((PyArrayObject*)array, NPY_ARRAY_ALIGNED | NPY_ARRAY_WRITEABLE);
+    return array;
+  }
+  %#endif
+}
 
 %typemap(in, fragment="fullMatrixConversion") const fullMatrix<double> &(PyObject *tmpObject = NULL, bool newMatrix = false){
   $1 = objToFullMatrixRO($input, newMatrix, tmpObject);
@@ -165,6 +313,19 @@
   if (newMatrix$argnum && $1) delete $1;
 }
 
+%typemap(in, fragment="fullMatrixConversionInt") const fullMatrix<int> &(PyObject *tmpObject = NULL, bool newMatrix = false){
+  $1 = objToFullMatrixIntRO($input, newMatrix, tmpObject);
+  if (!$1) {
+    PyErr_Format(PyExc_TypeError, "cannot convert argument %i to a fullMatrix<int>", $argnum);
+    SWIG_fail;
+  }
+}
+%typemap(freearg) const fullMatrix<int> &{
+  if (tmpObject$argnum) Py_DECREF(tmpObject$argnum);
+  if (newMatrix$argnum && $1) delete $1;
+}
+
+
 %typemap(in, fragment="fullMatrixConversion") fullMatrix<double> &(bool newMatrix = false){
   $1 = objToFullMatrixRW($input, newMatrix);
   if (!$1) {
@@ -176,6 +337,17 @@
   if (newMatrix$argnum && $1) delete $1;
 }
 
+%typemap(in, fragment="fullMatrixConversionInt") fullMatrix<int> &(bool newMatrix = false){
+  $1 = objToFullMatrixIntRW($input, newMatrix);
+  if (!$1) {
+    PyErr_Format(PyExc_TypeError, "cannot convert argument %i to a writable fullMatrix<int>", $argnum);
+    SWIG_fail;
+  }
+}
+%typemap(freearg) fullMatrix<int> &{
+  if (newMatrix$argnum && $1) delete $1;
+}
+
 %typemap(in, fragment="fullMatrixConversion") const fullMatrix<double> *(PyObject *tmpObject = NULL, bool newMatrix = false){
   $1 = objToFullMatrixRO($input, newMatrix, tmpObject);
   if (!$1) {
@@ -188,6 +360,18 @@
   if (newMatrix$argnum && $1) delete $1;
 }
 
+%typemap(in, fragment="fullMatrixConversionInt") const fullMatrix<int> *(PyObject *tmpObject = NULL, bool newMatrix = false){
+  $1 = objToFullMatrixIntRO($input, newMatrix, tmpObject);
+  if (!$1) {
+    PyErr_Format(PyExc_TypeError, "cannot convert argument %i to a fullMatrix<int>", $argnum);
+    SWIG_fail;
+  }
+}
+%typemap(freearg) const fullMatrix<int> *{
+  if (tmpObject$argnum) Py_DECREF(tmpObject$argnum);
+  if (newMatrix$argnum && $1) delete $1;
+}
+
 %typemap(in, fragment="fullMatrixConversion") fullMatrix<double> *(bool newMatrix = false){
   $1 = objToFullMatrixRW($input, newMatrix);
   if (!$1) {
@@ -199,6 +383,17 @@
   if (newMatrix$argnum && $1) delete $1;
 }
 
+%typemap(in, fragment="fullMatrixConversionInt") fullMatrix<int> *(bool newMatrix = false){
+  $1 = objToFullMatrixIntRW($input, newMatrix);
+  if (!$1) {
+    PyErr_Format(PyExc_TypeError, "cannot convert argument %i to a writable fullMatrix<int>", $argnum);
+    SWIG_fail;
+  }
+}
+%typemap(freearg) fullMatrix<int> *{
+  if (newMatrix$argnum && $1) delete $1;
+}
+
 #ifdef HAVE_NUMPY
 %typemap(out, fragment="fullMatrixConversion") fullMatrix<double> {
   $result = fullMatrix2PyArray($1);
@@ -215,5 +410,19 @@
 %typemap(out, fragment="fullMatrixConversion") fullVector<double>& {
   $result = fullVector2PyArrayProxy(*$1);
 }
+
+%typemap(out, fragment="fullMatrixConversionInt") fullMatrix<int> {
+  $result = fullMatrix2PyArray($1);
+}
+
+%typemap(out, fragment="fullMatrixConversionInt") const fullMatrix<int>& {
+  $result = fullMatrix2PyArrayConst(*$1);
+}
+
+%typemap(out, fragment="fullMatrixConversionInt") fullMatrix<int>& {
+  $result = fullMatrix2PyArrayProxy(*$1);
+}
+
 #endif
 
+