How to fix this?

when I run this code, it says that

Assertion fail : (th == &fe[1]->Vh->Th)
        line :4386, in file lgfem.cpp
Assertion fail : (th == &fe[1]->Vh->Th)
        line :4386, in file lgfem.cpp
 err code 6 ,  mpirank 0
terminate called after throwing an instance of 'std::length_error'
  what():  vector::_M_default_append
 try getConsole C:\Users\math\Desktop\Newton - copy.edp
 save log in :  'C:\Users\math\Desktop\Newton - copy.log'
wait enter ?

the complete code is

load "iovtk"
verbosity=0;
// Parameters
int nn = 15;
int nnPlus = 5;
real l = 1.;
real L = 15.;
real hSlope = 0.1;
real H = 6.;
real h = 0.5;

real beta = 0.01;

real eps = 9.81/303;
real nu = 1.;
real numu = nu/sqrt(0.09);
real nuep = pow(nu,1.5)/4.1;
real dt = 0.;

real Penalty = 1.e-6;

// Mesh
border b1(t=0, l){x=t; y=0;}
border b2(t=0, L-l){x=1.+t; y=-hSlope*t;}
border b3(t=-hSlope*(L-l), H){x=L; y=t;}
border b4(t=L, 0){x=t; y=H;}
border b5(t=H, h){x=0; y=t;}
border b6(t=h, 0){x=0; y=t;}

mesh Th=buildmesh(b1(nnPlus*nn*l) + b2(nn*sqrt((L-l)^2+(hSlope*(L-l))^2)) + b3(nn*(H + hSlope*(L-l))) + b4(nn*L) + b5(nn*(H-h)) + b6(nnPlus*nn*h));
plot(Th);

// Fespaces
fespace Vh2(Th, P1b);
Vh2 Ux, Uy,du1,du2,uite1,uite2;
Vh2 Vx, Vy,v1,v2;
Vh2 Upx, Upy=0;
fespace Vh(Th,P1);
Vh p=0, q,pp=p,pite,dp,vp;
Vh Tp, T=35,dT,Tite,vt;
Vh k=log(0.27), kp=k,dk,kite,vk;
Vh ep=log(0.002104), epp=ep,dep,epite,vep;
 
fespace V0h(Th,P0);
 V0h muT=1;
 V0h prodk, prode;
 Vh kappa=0.25e-4, stress;
// Macro
 macro grad(u) [dx(u), dy(u)] //
 macro Grad(u1,u2) [dx(u1),dy(u1), dx(u2),dy(u2)] //
 macro Gradt(u1,u2) [dx(u1),dx(u2), dy(u1),dy(u2)] //
 macro Div(u1,u2) (dx(u1) + dy(u2)) //
macro UgradV(u1,u2,v1,v2) [[u1,u2]'*[dx(v1),dy(v1)],[u1,u2]'*[dx(v2),dy(v2)]]//
macro ugradv(u1,u2,v) [u1,u2]'*grad(v)//
// Problem
real alpha = 0.;


// Initialization with stationary solution

plot([Upx, Upy],pp, value=true, coef=0.2, cmm="[Ux, Uy] ");
// Initialization
dt = 0.05; //probably too big
int nbiter = 1;
real coefdt = 0.25^(1./nbiter);
real coefcut = 0.25^(1./nbiter);
real cut = 0.01;
real tol = 0.5;
real coeftol = 0.5^(1./nbiter);

Tp = T - 10*((x<1)*(y<0.5) + (x>=1)*(y+0.1*(x-1)<0.5));
Upx=0+(x==0)*(y>=0)*(y<=h)*3;
// Convergence loop
real T0 = clock();
cout  << " - dt = " << dt << endl;
   alpha = 1/dt;

   // Time loop
   real t = 0.;
   for (int i = 0; i <= 500; i++){
      uite1=Upx;
      uite2=Upy;
      pite=pp;
      Tite=Tp;
      kite=kp;
      epite=epp;

      for (int n=0;n<=10;n++){
         solve increment ([du1,du2,dp,dT,dk,dep],[v1,v2,vp,vt,vk,vep])
         =int2d(Th)(alpha*([du1,du2]'*[v1,v2]-dT*vt-dk*vk-dep*vep)
         +0.5*(UgradV(du1,du2,uite1,uite2))'*[v1,v2]
         +0.5*(UgradV(uite1,uite2,du1,du2))'*[v1,v2]
         +0.5*ugradv(v1,v2,dp)
         +(0.09/2)*exp(2*kite-epite)*(2*dk-dep)*(Grad(uite1,uite2)+(Gradt(uite1,uite2)))'*Grad(v1,v2)
         +(0.09/2)*exp(2*kite-epite)*(Grad(du1,du2)+Gradt(du1,du2))'*Grad(v1,v2)
         +(9.81/606)*dT*v2-0.5*Div(du1,du2)*vp-0.5*vk*(ugradv(du1,du2,kite)+ugradv(uite1,uite2,dk))-0.5*exp(epite-kite)*(dep-dk)*vk+(0.09/2)*exp(2*kite-epite)*(2*dk-dep)*[dx(kite),dy(kite)]'*[dx(kite),dy(kite)]*vk
         +(0.09/2)*exp(2*kite-epite)*2*[dx(kite),dy(kite)]'*[dx(dk),dy(dk)]*vk
         -(0.09/2)*exp(2*kite-epite)*(2*dk-dep)*[dx(vk),dy(vk)]'*[dx(kite),dy(kite)]-(0.09/2)*exp(2*kite-epite)*[dx(vk),dy(vk)]'*[dx(dk),dy(dk)]
         +(0.09/2)*exp(kite-epite)*(4*dx(uite1)*dx(du1)+4*dy(uite2)*dy(du2)+2*(dy(uite1)+dx(uite2))*(dy(du1)+dx(du2)))*vk
         +(0.09/2)*exp(kite-epite)*(dk-dep)*(2*dx(uite1)*dx(uite1)+2*dy(uite2)*dy(uite2)+(dy(uite1)+dx(uite2))*(dy(uite1)+dx(uite2)))*vk
         -0.5*vep*(ugradv(du1,du2,epite)+ugradv(uite1,uite2,dep))
         -(1.92/2)*exp(-kite+epite)*(dep-dk)*vep
         +(0.09/2.6)*exp(2*kite-epite)*(2*dk-dep)*[dx(epite),dy(epite)]'*[dx(epite),dy(epite)]*vep
         +(0.09/2.6)*exp(2*kite-epite)*2*[dx(epite),dy(epite)]'*[dx(dep),dy(dep)]*vep
         -(0.09/2.6)*exp(2*kite-epite)*(2*dk-dep)*[dx(epite),dy(epite)]'*[dx(vep),dy(vep)]
         -(0.09/2.6)*exp(2*kite-epite)*[dx(dep),dy(dep)]'*[dx(vep),dy(vep)]
         +((1.44*0.09)/2)*exp(kite-epite)*(4*dx(uite1)*dx(du1)+4*dy(uite2)*dy(du2)+2*(dy(uite1)+dx(uite2))*(dy(du1)+dx(du2)))*vep
         +((1.44*0.09)/2)*exp(kite-epite)*(dk-dep)*(2*dx(uite1)*dx(uite1)+2*dy(uite2)*dy(uite2)+(dy(uite1)+dx(uite2))*(dy(uite1)+dx(uite2)))*vep
         -0.5*vt*(ugradv(du1,du2,kite)+ugradv(uite1,uite2,dk))
         -0.05*exp(2*kite-epite)*(2*dk-dep)*[dx(vt),dy(vt)]'*[dx(Tite),dy(Tite)]
         -0.05*exp(2*kite-epite)*[dx(vt),dy(vt)]'*[dx(dT),dy(dT)]
         )
         +int1d(Th,b1,b2,b5)(
            (0.09/2.6)*2.4952*exp(0.5*kite)*(0.5*dk)*vep
         )+int1d(Th,b1,b2)(
            0.41*((0.09)^0.25)*0.5*(exp(0.5*kite)*(0.5*dk)*uite1*v1+exp(0.5*kite)*du1*v1)
         )+int2d(Th)(
            alpha*([uite1-Upx,uite2-Upy]'*[v1,v2]-(Tite-Tp)*vt-(kite-kp)*vk-(epite-epp)*vep)
            +0.5*(UgradV(uite1,uite2,uite1,uite2))'*[v1,v2]+0.5*(UgradV(Upx,Upy,Upx,Upy))'*[v1,v2]
            +0.5*(grad(pite))'*[v1,v2]+0.5*(grad(pp))'*[v1,v2]
            +0.09*0.5*exp(2*kite-epite)*(Grad(uite1,uite2)+(Gradt(uite1,uite2)))'*Grad(v1,v2)
            +0.09*0.5*exp(2*kp-epp)*(Grad(Upy,Upy)+(Gradt(Upx,Upy)))'*Grad(v1,v2)
            +(9.81/606)*(Tite+Tp-70)*v2-0.5*(Div(uite1,uite2)+Div(Upx,Upy))*vp
            -0.5*(ugradv(uite1,uite2,kite))*vk-0.5*(ugradv(Upx,Upy,kp))*vk
            -0.5*exp(epite-kite)*vk-0.5*exp(epp-kp)*vk
            +(0.09/2)*vk*(exp(2*kite-epite)*grad(kite)'*grad(kite))
            +(0.09/2)*vk*(exp(2*kp-epp)*grad(kp)'*grad(kp))
            -(0.09/2)*(exp(2*kite-epite)*grad(kite)'*grad(vk))
            -(0.09/2)*(exp(2*kp-epp)*grad(kp)'*grad(vk))
            +(0.09/2)*exp(kite-epite)*vk*(2*dx(uite1)*dx(uite1)+2*dy(uite2)*dy(uite2)+(dy(uite1)+dx(uite2))*(dy(uite1)+dx(uite2)))
            +(0.09/2)*exp(kp-epp)*vk*(2*dx(Upx)*dx(Upx)+2*dy(Upy)*dy(Upy)+(dy(Upx)+dx(Upy))*(dy(Upx)+dx(Upy)))
            -0.5*(ugradv(uite1,uite2,epite))*vep-0.5*(ugradv(Upx,Upy,epp))*vep
            -1.92*0.5*exp(epite-kite)*vep-1.92*0.5*exp(epp-kp)*vep
            +(0.09/2.6)*vep*(exp(2*kite-epite)*grad(epite)'*grad(epite))+(0.09/2.6)*vep*(exp(2*kp-epp)*grad(epp)'*grad(epp))
            -(0.09/2.6)*(exp(2*kite-epite)*grad(epite)'*grad(vep))-(0.09/2.6)*(exp(2*kp-epp)*grad(epp)'*grad(vep))
            +(1.44*0.09/2)*exp(kite-epite)*vep*(2*dx(uite1)*dx(uite1)+2*dy(uite2)*dy(uite2)+(dy(uite1)+dx(uite2))*(dy(uite1)+dx(uite2)))
            +(1.44*0.09/2)*exp(kp-epp)*vep*(2*dx(Upx)*dx(Upx)+2*dy(Upy)*dy(Upy)+(dy(Upx)+dx(Upy))*(dy(Upx)+dx(Upy)))
            -0.5*(ugradv(uite1,uite2,Tite))*vt-0.5*(ugradv(Upx,Upy,Tp))*vt
            -0.05*(exp(2*kite-epite)*grad(Tite)'*grad(vt))-0.05*(exp(2*kp-epp)*grad(Tp)'*grad(vt))
         )+int1d(Th,b1,b2,b5)(
            (0.09/2.6)*2.4952*exp(0.5*kite)*vep+(0.09/2.6)*2.4952*exp(0.5*kp)*vep)
         +int1d(Th,b1,b2)(0.41*((0.09)^0.25)*exp(0.5*kite)*uite1*v1*0.5+0.41*((0.09)^0.25)*exp(0.5*kp)*Upx*v1*0.5)
         +on(b5,du1=0,du2=0)
         +on(b6,du1=0,du2=0,dep=0,dk=0,dT=0)+on(b2,dT=0)+on(b1,du2=0,dT=0);
         uite1=du1+uite1;
         uite2=du2+uite2;
         pite=dp+pite;
         Tite=dT+Tite;
         kite=dk+kite;
         epite=dep+epite;
      }
      Upx=uite1;
      Upy=uite2;
      pp=pite;
      Tp=Tite;
      kp=kite;
      epp=epite;
      plot([Upx, Upy],pp, coef=0.2, cmm=" [Ux, Uy] ", WindowIndex=1);

   }
cout << "Total Time = " << clock()-T0 << endl;

you could just ignore the long formula, I don’t know why this happen

I’m not sure if this helps but it did manage to get it to run past that point. If you look
at lgfem.cpp it looks like it is comparing the x and y pieces.

diff test1*
2d1
< load "medit"
32,33c31,32
< plot(Th,wait=true);
< medit("asdf",Th);
---
> plot(Th);
> 
38,40c37
< //Vh2 Upx, Upy=0;
< //Vh2 Upx, Upy;
< Vh2 Upx=0, Upy=0;
---
> Vh2 Upx, Upy=0;
42,43c39
< //Vh p=0, q,pp=p,pite,dp,vp;
< Vh p, q,pp=p,pite,dp,vp;
---
> Vh p=0, q,pp=p,pite,dp,vp;
64,66c60
< plot(Upx,cmm="Upx",wait=1);
< plot(Upy,cmm="Upy",wait=1);
< plot(pp,cmm="pp",wait=1);
---
> 
marchywka@happy:/home/documents/cpp/proj/freefem/play$ 

test1.edp (7.0 KB)

Upx is not defined, thank it is a bug, I will correct in next version
add Upx=0;