Use of PWMatrixCoefficient to define matrix coefficients per attribute

[mdavids-cfs]

Hi all, I have been trying to implement a PWMatrixCoefficient but I am having trouble with the interface (also couldn't find an example for this). The goal is to assign a MatrixConstantCoefficient per mesh attribute. In normal MFEM (C++) this is done by creating an Array Array<MatrixCoefficient *> coefs(0); and then assigning the matrix coefficients. I thought that PyMFEM can interpret normal lists of tuples as list-equivalents, i.e., I tried something like this:

sigma_all_coefs = []
sigma_attr = mfem.intArray(max_attr)
for ti,tensor in enumerate(sigma_all): 
    # just for testing
    temp = mfem.DenseMatrix(dim)
    temp.Assign(0.0)

    # add matrix coefficient to list
    sigma_all_coefs.append(mfem.MatrixConstantCoefficient(temp)  )
    sigma_attr[ti] = ti+1   
    
# Create PW Matrix Coefficient
sigmaCoef = mfem.PWMatrixCoefficient(dim, sigma_attr, sigma_all_coefs, False)

which fails when calling the mfem.PWMatrixCoefficient constructor. Also tried another approach:

# Create PW Matrix Coefficient
sigmaCoef = mfem.PWMatrixCoefficient(dim, False)

for ti, tensor in enumerate(sigma_all):
    # Set coefficient
    sigmaCoef.UpdateCoefficient( ti+1, mfem.MatrixConstantCoefficient( mfem.DenseMatrix(tensor) ) )

which fails when attempting to assemble a bilinear operator, that utilizes the PWMatrixCoefficient.

Any thoughts on how this could be achieved? Big thanks and warm regards, Mathias

[v-dobrev]

One option, as discussed in #57, may be to derive your own class from mfem.MatrixPyCoefficient as illustrated in https://github.com/mfem/PyMFEM/blob/master/examples/wrapper_example/matrix_coefficient.py.

Unfortunately, I don't know if one can create an equivalent of Array<MatrixCoefficient*> in python and I don't think sigma_all_coefs, as constructed in your first code snippet, has the right type. Also, in the second code snippet, when calling sigmaCoef.UpdateCoefficient the second argument in C++ is MatrixCoefficient & and I suspect that on the python side using mfem.MatrixConstantCoefficient is not correct but I don't know why -- I'm not familiar with how SWIG translates things.

[sshiraiwa]

@mdavids-cfs, @v-dobrev This is going to be address in https://github.com/mfem/PyMFEM/tree/coefficient-arrray-dev, where MatrixCoefficientPtrArray is available. The second option may work in the current master, if you keep the argument objects ( mfem.MatrixConstantCoefficient and mfem.DenseMatrix) from being gc-ed.

[mdavids-cfs]

Hi @v-dobrev and @sshiraiwa I should've responded earlier; I did in fact get it to working by storing MatrixConstantCoefficient and DenseMatrix in lists, so they don't get overwritten, and then using the second approach. This worked pretty well! Also, big thanks for putting together this python interface to MFEM, it's really helpful!