Segmentation fault when assembling linear form with MatrixVectorProductCoefficient

[benzwick]

I am trying to use MatrixVectorProductCoefficient in a DomainLFGradIntegrator as follows but I get a segmentation fault:

With Numba:

import mfem.ser as mfem
import numpy as np
from numba import cfunc, carray

mesh = mfem.Mesh(10, 10, 10, "HEXAHEDRON")
sdim = mesh.SpaceDimension()
ndim = mesh.Dimension()

fec = mfem.H1_FECollection(1, mesh.Dimension())
fespace = mfem.FiniteElementSpace(mesh, fec)

k = 2.0
C = mfem.MatrixConstantCoefficient([[ k, 0., 0.],
                                    [0.,  k, 0.],
                                    [0., 0.,  k]])

@cfunc(mfem.vector_sig)
def v_func(x, out, sdim, vdim):
    out_array = carray(out, (vdim, ))
    out_array[0] = x[0]
    out_array[1] = x[1]
    out_array[2] = x[2]

v_coef = mfem.VectorNumbaFunction(v_func, sdim, ndim).GenerateCoefficient()
f = mfem.MatrixVectorProductCoefficient(C, v_coef)

b = mfem.LinearForm(fespace)
b.AddDomainIntegrator(mfem.DomainLFGradIntegrator(f))
b.Assemble()

Without Numba:

import mfem.ser as mfem
import numpy as np

mesh = mfem.Mesh(10, 10, 10, "HEXAHEDRON")
ndim = mesh.Dimension()

fec = mfem.H1_FECollection(1, mesh.Dimension())
fespace = mfem.FiniteElementSpace(mesh, fec)

k = 2.0
C = mfem.MatrixConstantCoefficient([[ k, 0., 0.],
                                    [0.,  k, 0.],
                                    [0., 0.,  k]])

class v_func(mfem.VectorPyCoefficient):
    def __init__(self):
       mfem.VectorPyCoefficient.__init__(self, ndim)
    def EvalValue(self, x):
       return (x[1], x[2], x[0])

v_coef = v_func()
f = mfem.MatrixVectorProductCoefficient(C, v_coef)

b = mfem.LinearForm(fespace)
b.AddDomainIntegrator(mfem.DomainLFGradIntegrator(f))
b.Assemble()

Does anyone know what is going wrong here?

EDIT: simplified the example and added non-Numba version.

[sshiraiwa]

@benzwick Thank you for reporting this. Looking into this, it seems this odd error comes from mfem itself, not wrapper.

The error is happening specifically this line, https://github.com/mfem/mfem/blob/master/fem/coefficient.cpp#L628

In mfem4.2, which PyMFEM is built with, this SetSize does not exist and the vector size remains zero.

If you are building mfem using setup.py, you can add SetSize to mfem and build it. To do this, please go to PyMFEM/external/mfem. This is a folder mfem is downloaded. You can edit it and in the PyMFEM base directory, run this command.

python setup.py install --ext-only (--with-parallel if you need this too)

Probably, this was a bug fixed between mfem-4.2 and 4.3. Therefore another option is to use mfem-4.3-dev of PyMFEM.
You can git clone mfem-4.3-dev branch of PyMFEM and use

python3 setup.py install --mfem-branch mfem-4.3-dev

You can find more discussion about this branch in #87

[benzwick]

I installed the mfem-4.3-dev branch but now I get the following error when creating the MatrixConstantCoefficient:

      9 
     10 k = 2.0
---> 11 C = mfem.MatrixConstantCoefficient([[ k, 0., 0.],
     12                                     [0.,  k, 0.],
     13                                     [0., 0.,  k]])

~/.local/lib/python3.8/site-packages/mfem/_ser/coefficient.py in __init__(self, m)
    789 
    790         if can_np_array:
--> 791            v = mfem._ser.vector.Vector(np.transpose(value).flatten())
    792            m = mfem._ser.densemat.DenseMatrix(v.GetData(), value.shape[0], value.shape[1])
    793            self._value = (v,m)

~/.local/lib/python3.8/site-packages/mfem/_ser/vector.py in __init__(self, *args)
    469 
    470 
--> 471         _vector.Vector_swiginit(self, _vector.new_Vector(*args))
    472 
    473         if keep_link:

TypeError: Wrong number or type of arguments for overloaded function 'new_Vector'.
  Possible C/C++ prototypes are:
    mfem::Vector::Vector()
    mfem::Vector::Vector(mfem::Vector const &)
    mfem::Vector::Vector(int)
    mfem::Vector::Vector(double *,int)
    mfem::Vector::Vector(int,mfem::MemoryType)
    mfem::Vector::Vector(int,mfem::MemoryType,mfem::MemoryType)
    mfem::Vector::Vector(mfem::Vector const &,int,int)

A similar error occurs if I try to create a vector:

In [1]: v = mfem.Vector([1.,2.,3.])                                                                 
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-2-5d39f275047f> in <module>
----> 1 v = mfem.Vector([1.,2.,3.])

~/.local/lib/python3.8/site-packages/mfem/_ser/vector.py in __init__(self, *args)
    469 
    470 
--> 471         _vector.Vector_swiginit(self, _vector.new_Vector(*args))
    472 
    473         if keep_link:

TypeError: Wrong number or type of arguments for overloaded function 'new_Vector'.
  Possible C/C++ prototypes are:
    mfem::Vector::Vector()
    mfem::Vector::Vector(mfem::Vector const &)
    mfem::Vector::Vector(int)
    mfem::Vector::Vector(double *,int)
    mfem::Vector::Vector(int,mfem::MemoryType)
    mfem::Vector::Vector(int,mfem::MemoryType,mfem::MemoryType)
    mfem::Vector::Vector(mfem::Vector const &,int,int)

[benzwick]

It appears that the bug in MatrixVectorProductCoefficient was reported in mfem/mfem#1976, fixed in mfem/mfem#1977 and merged into the master branch that will become 4.3.0.

Do you know if there are any plans to release MFEM 4.2.1 etc. with fixes for bugs like this?

[sshiraiwa]

MFEM 4.3 seems coming soon. We should better just fix the issue you mentioned above. Let me look into it.

[benzwick]

OK, MatrixVectorProductCoefficient works with 4.2 if I make the change to coefficient as you described above so shall I close this issue and open a new one for the problem with vector in 4.3 or just leave this one open for now?

[sshiraiwa]

The issue in mfem-4.3 is due to the recent change of mfem header, where _data (and any other variable starting
with underscore) is renamed.

[sshiraiwa]

Closing this for now.
Fix was made in current mfem-4.3-dev. Please re-open if you have further problems.