Creating Custom Matrices and Vectors from Example 25 Code Results in Memory Leaks

[cyasing]

In Electromagnetics, I tried Example 25 and it runs perfectly fine (with GLVis plots). My goal is to be able to create new matrices from the default ax = b system for further processing. However, no matter what I do, I seem to encounter some form of exception due to memory not being accessible / not being freed in time.

For instance, this is a map of DOFs from global to local numbering that I am trying to write just after you recover the solution in example 25. I am trying for problem type iprob = 4:

//    FiniteElementSpace* fespace_x = x.FESpace();
//    int num_dofs_x = fespace_x->GetVSize();
//    Array<int> dof_map_x(num_dofs_x); // Mapping from global DOF indices to local node indices
//
//    for (int i = 0; i < fespace_x->GetNE(); i++)
//    {
//        Array<int> dofs;
//        fespace_x->GetElementDofs(i, dofs); // Get local DOFs for this element
//        for (int j = 0; j < dofs.Size(); j++)
//        {
//            dof_map_x[dofs[j]] = fespace_x->DofToVDof(i, j); // Map local DOFs to global DOFs
//        }
//    }

I tried deleting the fespace_x in any order with respect to the original list of delete, but doesn't work. I understand since we are not allocating memory using new, we do not really need delete, but I still get an exception even without specifying delete fespace_x (right at return 0;).

Similar things happen when I try to create a separate matrix for just the time derivative term (-\omega^2 * epsilon E):

    VectorFEMassIntegrator* m = new VectorFEMassIntegrator(restr_omeg);
    SesquilinearForm massTerm(fespace, conv);
    massTerm.AddDomainIntegrator(m, NULL);
    massTerm.Assemble(0);
    SparseMatrix& M_real = massTerm.real().SpMat(); // Real part of the mass matrix

Please help me understand how such matrices should be constructed.

[psocratis]

Hi @cyasing

The main issue here is that dofs[j] might return a negative index depending on the orientation of the ND dofs. To decode it you can use FiniteElementSpace::DecodeDof(dofs[j]) to get the non-negative index. I suspect that that's why you are getting the memory exception. In addition, you don't need the DofToVdof since this is an ND space with vdim=1 (vdim indicates the number of copies of a scalar finite element space, such as H1^dim). If you want to get the element dof to global dof map this would simply be j->FiniteElementSpace::DecodeDof(dofs[j])

I don't see an issue with the construction of the mass matrix.
If you keep running into issues, feel free to post the whole driver that reproduces the error. I also advise compiling MFEM in debug mode. This can give more information on potential issues.

I hope this helps.

Socratis

[cyasing]

Thanks for your response! The DOF mapping works now, but I still get an error for the mass matrix. The exception happens in bilinear.cpp at line 1186 for (k=0; k < domain_integs.Size(); k++) { delete domain_integs[k]; }. A few points I have doubts on:

1. Where should I place the code I have mentioned above for the mass matrix? Currently it's below the line where the linear system is formed in the example: a.FormLinearSystem(ess_tdof_list, x, b, A, X, B); altogether.
2. Should the assembly of this matrix be kept before the assembly of the whole linear system? Or should it be a part of assembly of a in some way, so that I don't create more than one full system?

[psocratis]

@cyasing,

The issue is that you are trying to delete the integrator. You shouldn't do this because the ownership of the integrator is passed to the bilinearform, i.e, the bilinearform will delete it.

It doesn't matter where you place it. It's independent of the rest of the code.

Best,

Socratis

[cyasing]

But here's the problem, I am not deleting anything other than the default deletes in the example. Which is why I am confused as to where the memory issue is originating from. This is all I delete:

    delete pml;
    delete fec;
    delete fespace;
    delete mesh;

[psocratis]

are you using the pointer m anywhere else?

[cyasing]

Ah shucks, you're right! I was using m to create the initial a matrix, and then reused it to create the mass matrix. But now I undid the pointer m declaration and separated the integrator for both matrices, and it works. I thought that since its the same term I could reuse the pointer instead of trying to create another new integrator. But why is it not the case?

[psocratis]

That's because the pointer to the integrator was owned by two different bilinear forms. When they went out of scope they both tried to delete the integrator.

Glad it works now.

[cyasing]

That makes sense, thanks!