# How to fix this?

**URL:** <https://community.freefem.org/t/how-to-fix-this/2327>\
**Category:** General Discussion\
**Created:** [March 12, 2023, 7:17pm UTC](https://community.freefem.org/t/how-to-fix-this/2327 "2023-03-12T19:17:12Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![cfdyuewu](https://avatars.discourse-cdn.com/v4/letter/c/5f8ce5/32.png) [@cfdyuewu](https://community.freefem.org/u/cfdyuewu)\
**Post date:** [March 12, 2023, 7:17pm UTC](https://community.freefem.org/t/how-to-fix-this/2327/1 "2023-03-12T19:17:12Z")

</div>

when I run this code, it says that

```auto
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

```auto
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

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [March 12, 2023, 8:29pm UTC](https://community.freefem.org/t/how-to-fix-this/2327/2 "2023-03-12T20:29:45Z")

</div>

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.

```auto
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](https://community.freefem.org/uploads/short-url/vShUSdFcCRXrSps5eCQe3Qp9raR.edp) (7.0 KB)

---

<div class="post-metadata">

**Author:** ![frederichecht](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/frederichecht/32/15_2.png) [@frederichecht](https://community.freefem.org/u/frederichecht)\
**Post date:** [March 13, 2023, 1:56pm UTC](https://community.freefem.org/t/how-to-fix-this/2327/3 "2023-03-13T13:56:45Z")

</div>

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