VectorPyCoefficent invalid size causing crash

[ddkn]

Hello everyone,

I was too quick to praise success in Issue 184 with

This makes getting stress trivial, now that I understand how this works better. Now on to more fun!

I have now created a stress vector coefficient as such

class StressVectorCoefficient(mfem.VectorPyCoefficient):
    def __init__(
            self, 
            lambda_coef: mfem.PWConstCoefficient, 
            mu_coef: mfem.PWConstCoefficient,
            dim: int):
        super(StressVectorCoefficient, self).__init__(dim)
        self.dim = dim
        self.lam = lambda_coef   # coefficient
        self.mu = mu_coef       # coefficient
        self.u = None   # displacement GridFunction
        self.grad = mfem.DenseMatrix()
        self.sigma = mfem.DenseMatrix()

    def SetDisplacement(self, u: mfem.GridFunction):
        self.u = u

    def Eval(
            self, v: mfem.Vector, 
            T: mfem.ElementTransformation, 
            ip: mfem.IntegrationPoint):
        L = self.lam.Eval(T, ip)
        M = self.mu.Eval(T, ip)
        
        self.u.GetVectorGradient(T, self.grad)
        self.grad.Symmetrize()
        
        print('g', self.grad.Size())
        # sigma = lambda * div_u * I = lambda * trace(grad(u)) * I
        self.sigma.Diag(L * self.grad.Trace(), self.grad.Size())
        print('a')
        # sigma += 2 * mu * grad(u)
        self.sigma.Add(2 * M, self.grad)
        print('b')
        
        size = self.sigma.Size()
        print('s', size)
        # Set the principle axes xx, yy, zz, xy, xz, yz
        if size == 2:
            v[0] = self.sigma[0, 0]
            v[1] = self.sigma[1, 1]
            v[2] = self.sigma[0, 1]
        elif size == 3:
            v[0] = self.sigma[0, 0]
            v[1] = self.sigma[1, 1]
            v[2] = self.sigma[2, 2]
            v[3] = self.sigma[0, 1]
            v[4] = self.sigma[0, 2]
            v[5] = self.sigma[1, 2]

        return self.sigma 

stress_coef = StressVectorCoefficient(lamb_coef, mu_coef, dim)
stress_coef.SetDisplacement(x)

vector_space = mfem.FiniteElementSpace(mesh, fec, dim)
stress = mfem.GridFunction(vector_space)
stress.ProjectCoefficient(stress_coef) 

The program only runs for a few iterations printing out

g 3
a
b
s 3
g 3
a
b
s 3
g 3
a
b
s 3
g 3
a
b

Meaning at some point it fails to get the self.sigma.Size(), the output from the console is this

...
   Iteration : 499  (B r, r) = 0.0244504
   Iteration : 500  (B r, r) = 0.0174414
PCG: Number of iterations: 500
Average reduction factor = 0.986248
PCG: No convergence!
free(): invalid next size (fast)

It can also read realloc(): invalid next size. Has anyone else dealt with this issue? For context in my other issue I calculate strain similarly but only do this and it works fine :/, the same model, everything

...
    def Eval(self, v, T, ip):
        self.u.GetVectorGradient(T, self.grad)
        self.grad.Symmetrize()
        
        size = self.grad.Size()
        # Set the principle axes xx, yy, zz, xy, xz, yz
        if size == 2:
            v[0] = self.grad[0, 0]
            v[1] = self.grad[1, 1]
            v[2] = self.grad[0, 1]
        elif size == 3:
            v[0] = self.grad[0, 0]
            v[1] = self.grad[1, 1]
            v[2] = self.grad[2, 2]
            v[3] = self.grad[0, 1]
            v[4] = self.grad[0, 2]
            v[5] = self.grad[1, 2]

[ddkn]

This now seems to be a problem with my original StrainVectorCoefficient which is leading me to think something is wrong with the mesh creating nodes that are invalid? I have no idea how to troubleshoot this.

I generate meshes on gmsh, Here is my script:
output.geo.txt

I have a python script that changes a few lines to change where the boundary is, which seems to work fine and renders in GLVis.

[ddkn]

Hm, I am noticing I need to add a import time; time.sleep(3) for this to work... now for just my StrainVectorCoefficient code. Is there a way to enforce a calculation is done before continuing.

Error

...
   Iteration : 135  (B r, r) = 0.000272998
   Iteration : 136  (B r, r) = 0.000214033
   Iteration : 137  (B r, r) = 0.000154586
   Iteration : 138  (B r, r) = 0.000109668
   Iteration : 139  (B r, r) = 7.95231e-05
Average reduction factor = 0.935443
Setting Nodal FE Space
Assigning reference nodes
data
Strain
[5.851080012110099e-06, -3.056506058355944e-06, -4.061944796849572e-06, 1.0322523091085359e-05]
corrupted double-linked list
Aborted
$ 

But I also get intermittent errors such as,

   Iteration : 136  (B r, r) = 0.000214033
   Iteration : 137  (B r, r) = 0.000154586
   Iteration : 138  (B r, r) = 0.000109668
   Iteration : 139  (B r, r) = 7.95231e-05
Average reduction factor = 0.935443
Setting Nodal FE Space
Assigning reference nodes
munmap_chunk(): invalid pointer
Aborted

Strain

class StrainVectorCoefficient(mfem.VectorPyCoefficient):                                                                                                                             
    def __init__(                                                                                                                                                                    
            self,                                                                                                                                                                    
            dim: int):                                                                                                                                                               
        super(StrainVectorCoefficient, self).__init__(dim)                                                                                                                           
        self.dim = dim                                                                                                                                                               
        self.u = None   # displacement GridFunction                                                                                                                                  
        self.grad = mfem.DenseMatrix()                                                                                                                                               
       
   def SetDisplacement(self, u: mfem.GridFunction):
        self.u = u

    def Eval(
            self, v: mfem.Vector, 
            T: mfem.ElementTransformation, 
            ip: mfem.IntegrationPoint) -> None:
        self.u.GetVectorGradient(T, self.grad)
        self.grad.Symmetrize()
        
        size = self.grad.Size()
        # Set the principle strains xx, yy, zz, xy, xz, yz
        if size == 2:
            v[0] = self.grad[0, 0]
            v[1] = self.grad[1, 1]
            v[2] = self.grad[0, 1]
        elif size == 3:
            v[0] = self.grad[0, 0]
            v[1] = self.grad[1, 1]
            v[2] = self.grad[2, 2]
            v[3] = self.grad[0, 1]
            v[4] = self.grad[0, 2]
            v[5] = self.grad[1, 2]

        return v


def get_xyz(c, mesh, gf):
    _, elem_ids, ips = mesh.FindPoints([c])
    c = [gf.GetValue(elem_ids[0], ips[0], i) for i in range(4)]
    return c

dim = 3
strain_coef = StrainVectorCoefficient(dim)
strain_coef.SetDisplacement(x)

vector_space = mfem.FiniteElementSpace(mesh, fec, dim)
strain = mfem.GridFunction(vector_space)
strain.ProjectCoefficient(strain_coef)    

print("data")
# WARN: if sleep is lower than 3 it crashes prematurely
sleep(5)
print("Strain", get_xyz((10, 0, 5), mesh, strain))

# NOTE: Complains about destructors missing
del strain
del strain_coef

The original issue is still present with StressVectorCoefficient. I can work with strain, but I am getting very intermittent errors on even if I can calculate the strain. So it isn't reliable.

[ddkn]

Hm, Opening in Paraview and moving shows some tetragonal elements missing as I rotate the solution... is it possible that mfem is dropping elements and unable to calculate the mesh?

A more basic version of the mesh, as two simple cylinders.

If I stop moving it, it smooths it out which makes me think interpolation is happening to fill in the missing regions?

So rebuilding the geo file to ensure it is coherent has got me back to calculating strain, via VectorPyCoefficient but it works only sometimes. I am back to the issue with StressVectorCoefficient generating the error:

free(): invalid next size (fast)

[sshiraiwa]

Hi @ddkn,.
This visualization looks strange to me, although I don't use paraview. Also, I am not sure about the sleep thing. Can you put a short program which demonstrates your problem so that I can take a close look? Thank you.

[ddkn]

Hi @sshiraiwa, sorry for getting back late, grad school has been a bit chaotic. Here is an example with a script, and GMSH geofile and msh, example.zip. The StressVectorCoefficient still continues to segfault,

One issue I discovered, is if I try to measure a spot that does not exist, the system crashes. I would likely need to implement some interpolation. The strain vector should work now. I did notice that small changes in the mesh generation can wildly change the element number.

As of now I would say being able to render via glvis -m test.mesh -g stress.gf would be considered a wild success, haha.

Thanks again for taking your time to help me out!