Error during matrix-vector multiplication with ffddm

General Discussion
1

Using ffddm, following the example Maxwell-3d-simple.edp in the examples/ffddm folder, I am trying to see what each operation does by printing out the results.
So here I am multiplying the partition of unity matrix pr#Dih[mpiRank(pr#commddm)] with a vector. I have defined x1 to be:

complex[int] x1(rhs.n);
x1 = 1;

The error lies in this line:

x1 = MDih[mpiRank(Mcommddm)] * x1;

with the error:
Assertion fail : ((!this->N() || !this->M()) || !Ax.x.N() || !this->N() || &Ax.x[0] != &this->operator)
line :560, in file ./…/femlib/RNM.hpp
err code 6
I have checked that the output of MDih[mpiRank(Mcommddm)] .n and x1.n are the same so it shouldn’t be a problem of unmatching dimension. But I don’t really understand what the error is telling me.
Below is my code for reference:

load "iovtk"
load "Element_Mixte3d"

macro dimension 3// EOM

include "ffddm.idp"

load "msh3"

macro def(i)[i, i#y, i#z]// EOM             // vector field definition
macro init(i)[i, i, i]// EOM                // vector field initialization
macro defPart(u)u// EOM                     // partition of unity definition
macro initPart(u)u// EOM                    // partition of unity initialization
macro Curl(ux, uy, uz)[dy(uz)-dz(uy), dz(ux)-dx(uz), dx(uy)-dy(ux)]// EOM
macro CrossN(ux, uy, uz)[uy*N.z-uz*N.y, uz*N.x-ux*N.z, ux*N.y-uy*N.x]// EOM
macro Curlabs(ux, uy, uz)[abs(dy(uz)-dz(uy)), abs(dz(ux)-dx(uz)), abs(dx(uy)-dy(ux))]// EOM
func Pk = Edge03d;
func PkPart = Edge03ds0;

int Dirichlet = 1, Robin = 2;
int ra = getARGV("-ra", 0);
int ng = getARGV("-global", 30);
int bc = getARGV("-bc", 1);
real nk = getARGV("-k", 10.);
real mul = getARGV("-mul", 40.);

int[int] labs(1);
if(bc == 1){labs = [Robin, Robin, Robin, Robin, Robin, Robin];}
else{labs = [Dirichlet, Dirichlet, Dirichlet, Dirichlet, Dirichlet, Dirichlet];}

// labs: list of surface labels, default=[1, 2, 3, 4, 5, 6]
mesh3 Th = cube(ng, ng, ng, [x,y,z], label = labs);

func f = exp(-((x-0.5)^2+(y-0.5)^2+(z-0.5)^2)*mul);
int mysplit = 2;

mesh3 Thc = cube(ng/mysplit, ng/mysplit, ng/mysplit, [x,y,z], label = labs);
// mesh3 Thc = cube(7, 7, 7, [x,y,z], label = labs);
// func k = 6*pi*nk;
func k = nk;
real epsilonprob = k^2;//0 k k^2;                  // epsilon of the problem
real be = getARGV("-betaEprec",2.); // beta for epsilonEprec
real epsilonEprec = k^be;//k^be;

//macro Mdefplot(u)sqrt(real(u)^2+real(u#y)^2+real(u#z)^2)//

vsym = 2; // symmetric, not hermitian
vtgv = 1.e+30;
vtgvelim = 1.e+30;

ffddmbuildDmesh(M,Th,mpiCommWorld)
ffddmbuildDfespaceEdge(M,M,complex,def,init,Pk,defPart,initPart,PkPart)

macro Varf(varfName, meshName, PhName)
    varf varfName([Ex,Ey,Ez],[vx,vy,vz]) =
  int3d(meshName)(Curl(vx,vy,vz)'*Curl(Ex,Ey,Ez))
                + int3d(meshName)(-(k^2-1i*epsilonprob)*[vx,vy,vz]'*[Ex,Ey,Ez])
                + int2d(meshName,Robin)(1i*k*CrossN(vx,vy,vz)'*CrossN(Ex,Ey,Ez))
                + on(Dirichlet,Ex=0,Ey=0,Ez=0);    // EOM

// for the preconditioner
macro VarfEprec(varfName, meshName, PhName)
    varf varfName([Ex,Ey,Ez],[vx,vy,vz]) =
  int3d(meshName)(Curl(vx,vy,vz)'*Curl(Ex,Ey,Ez))
                + int3d(meshName)(-(k^2-1i*epsilonEprec)*[vx,vy,vz]'*[Ex,Ey,Ez])
                + int2d(meshName,Robin)(1i*k*CrossN(vx,vy,vz)'*CrossN(Ex,Ey,Ez))
                + on(Dirichlet,Ex=0,Ey=0,Ez=0);    // EOM
                
// for the preconditioner
macro VarfOpt(varfName, meshName, PhName)

    varf varfName([Ex,Ey,Ez],[vx,vy,vz]) =
  int3d(meshName)(Curl(vx,vy,vz)'*Curl(Ex,Ey,Ez))
                + int3d(meshName)(-(k^2-1i*epsilonEprec)*[vx,vy,vz]'*[Ex,Ey,Ez])
                + int2d(meshName,Robin)(1i*k*CrossN(vx,vy,vz)'*CrossN(Ex,Ey,Ez))
                + int2d(meshName,10)(1i*k*CrossN(vx,vy,vz)'*CrossN(Ex,Ey,Ez))
                + on(Dirichlet,Ex=0,Ey=0,Ez=0);    // EOM

macro Varfrhs(varfName, meshName, PhName)
    varf varfName([Ex,Ey,Ez],[vx,vy,vz]) =
        -int3d(meshName)([vx,vy,vz]'*[0,0,f])
       + on(Dirichlet,Ex=0,Ey=0,Ez=0);  // EOM


ffddmsetup(M,M,Varf,Varfprec)
if (ra == 1) {ffddmsetup(M, M, Varf, VarfOpt)}

complex[int] rhs(1);

ffddmbuildrhs(M,Varfrhs,rhs)

complex[int] x0(rhs.n);
x0 = 0;

MVhi<complex> def(u);

u[] = MfGMRES(x0, rhs, 1.e-6, 800, "right");

complex[int] x1(rhs.n);
x1 = 1;

cout << "Process " << mpiRank(Mcommddm) << "'s Di has " << MDih[mpiRank(Mcommddm)].n << " dofs." << endl;
x1 = MDih[mpiRank(Mcommddm)] * x1;
2

Unfortunately you cannot use the same vector as input and output for such operations in FreeFEM. You need to define a second vector as output :

complex[int] x1(rhs.n), x2(rhs.n);
x1 = 1;

x2 = MDih[mpiRank(Mcommddm)] * x1;

or directly

complex[int] x1(rhs.n);
x1 = 1;

complex[int] x2 = MDih[mpiRank(Mcommddm)] * x1;
3

Thanks!
It works now.
By the way,
I found that the elementwise vector multiplication x1 = x1 .* MDk[mpiRank(Mcommddm)]; works too. It seems a bit mysterious to me why the matrix-vector product won’t work.
But thanks for the reply! It solved my problem.