# How to set the initial value on the boundary?

I’m working on the example “large fluid problem” by using fem and newton iteration method, I want to set the value of Ux on b6 be 3, because newton method need to solve the incremental value, so I can’t directly assign the dirichlet boundary condition with respect to Ux,Uy,here is the code

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=-hSlopet;}
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(nnPlusnnl) + b2(nnsqrt((L-l)^2+(hSlope(L-l))^2)) + b3(nn*(H + hSlope*(L-l))) + b4(nnL) + b5(nn(H-h)) + b6(nnPlusnnh));
plot(Th);

// Fespaces
fespace Vh2(Th, P1b);
Vh2 Ux, Uy,dux,duy;
Vh2 Vx, Vy;
for (int i = 0; i < Th.nv; ++i)
{
if (Th(i).onboundary && Th(i).boundary == 6)
Ux[i] = 3.0;
}
Ux=0+on(b6,Ux=3);
Vh2 Upx=Ux, Upy=Uy;
fespace Vh(Th,P1);
Vh p=0, q,pp;
Vh Tp, T=35,dT;
Vh k=log(0.27), kp=k,dk;
Vh ep=log(0.002104), epp=ep,dep;

fespace V0h(Th,P0);
V0h muT=1;
V0h prodk, prode;
Vh kappa=0.25e-4, stress;
// Macro
macro Div(U) (dx(U#x) + dy(U#y)) //
// Functions
func g = (x) * (1-x) * 4;
// Problem
real alpha = 0.;

// Initialization with stationary solution

plot([Ux, Uy], p, value=true, coef=0.2, cmm=“[Ux, Uy] - p”);
// 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);

T = T - 10*((x<1)(y<0.5) + (x>=1)(y+0.1*(x-1)<0.5));

// Convergence loop
real T0 = clock();
for (int iter = 1; iter <= nbiter; iter++){
cout << "Iteration " << iter << " - dt = " << dt << endl;
alpha = 1/dt;

// Time loop
real t = 0.;
for (int i = 0; i <= 500; i++){
u_ite1=Upx;
u_ite2=Upy;
p_ite=pp;
T_ite=Tp;
k_ite=kp;
ep_ite=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.09/2)exp(2*k_ite-ep_ite)*2*[dx(k_ite),dy(k_ite)]'*[dx(dk),dy(dk)]*vk
-(0.09/2)exp(2*k_ite-ep_ite)*(2*dk-dep)*[dx(vk),dy(vk)]'*[dx(k_ite),dy(k_ite)]-(0.09/2)exp(2*k_ite-ep_ite)*[dx(vk),dy(vk)]'*[dx(dk),dy(dk)]
+(0.09/2)*exp(k_ite-ep_ite)*(4*dx(u_ite1)*dx(du1)+4*dy(u_ite2)*dy(du2)+2*(dy(u_ite1)+dx(u_ite2))*(dy(du1)+dx(du2)))*vk
+(0.09/2)*exp(k_ite-ep_ite)*(dk-dep)*(2*dx(u_ite1)*dx(u_ite1)+2*dy(u_ite2)*dy(u_ite2)+(dy(u_ite1)+dx(u_ite2))*(dy(u_ite1)+dx(u_ite2)))*vk
-(1.92/2)*exp(-k_ite+ep_ite)*(dep-dk)*vep
+(0.09/2.6)*exp(2*k_ite-ep_ite)*(2*dk-dep)*[dx(ep_ite),dy(ep_ite)]'*[dx(ep_ite),dy(ep_ite)]*vep
+(0.09/2.6)*exp(2*k_ite-ep_ite)*2*[dx(ep_ite),dy(ep_ite)]'*[dx(dep),dy(dep)]*vep
-(0.09/2.6)*exp(2*k_ite-ep_ite)*(2*dk-dep)*[dx(ep_ite),dy(ep_ite)]'*[dx(vep),dy(vep)]
-(0.09/2.6)*exp(2*k_ite-ep_ite)*[dx(dep),dy(dep)]'*[dx(vep),dy(vep)]
+((1.44*0.09)/2)*exp(k_ite-ep_ite)*(4*dx(u_ite1)*dx(du1)+4*dy(u_ite2)*dy(du2)+2*(dy(u_ite1)+dx(u_ite2))*(dy(du1)+dx(du2)))*vep
+((1.44*0.09)/2)*exp(k_ite-ep_ite)*(dk-dep)*(2*dx(u_ite1)*dx(u_ite1)+2*dy(u_ite2)*dy(u_ite2)+(dy(u_ite1)+dx(u_ite2))*(dy(u_ite1)+dx(u_ite2)))*vep
-0.05*exp(2*k_ite-ep_ite)*(2*dk-dep)*[dx(vt),dy(vt)]'*[dx(T_ite),dy(T_ite)
-0.05*exp(2*k_ite-ep_ite)*[dx(vt),dy(vt)]'*[dx(dT),dy(dT)
)
+int1d(Th,b1,b2,b5)(
(0.09/2.6)*2.4952*exp(0.5*k_ite)*(0.5*dk)*vep
)+int1d(Th,b1,b2)(
0.41*((0.09).^0.25)*0.5*(exp(0.5*k_ite)*(0.5*dk)*u_ite1*v1+exp(0.5*k_ite)*du1*v1)
)+int2d(Th)(
alpha*([u_ite1-Upx,u_ite2-Upy]'*[v1,v2]-(T_ite-Tp)*vt-(k_ite-kp)*vk-(ep_ite-epp)*vep)
+(9.81/606)*(T_ite+Tp-70)*v2-0.5*Div([u_ite1+Upx,u_ite2+Upy])*vp
-0.5*exp(ep_ite-k_ite)*vk-0.5*exp(epp-kp)*vk
+(0.09/2)*exp(k_ite-ep_ite)*vk*(2*dx(u_ite1)*dx(u_ite1)+2*dy(u_ite2)*dy(u_ite2)+(dy(u_ite1)+dx(u_ite2))*(dy(u_ite1)+dx(u_ite2)))
+(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)))
-1.92*0.5*exp(-k_ite+ep_ite)*vep--1.92*0.5*exp(-kp+epp)*vep
+(1.44*0.09/2)*exp(k_ite-ep_ite)*vep*(2*dx(u_ite1)*dx(u_ite1)+2*dy(u_ite2)*dy(u_ite2)+(dy(u_ite1)+dx(u_ite2))*(dy(u_ite1)+dx(u_ite2)))
+(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)))
)+int1d(Th,b1,b2,b5)(
(0.09/2.6)*2.4952*exp(0.5*k_ite)*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*k_ite)*u_ite1*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, du2=-du1*N.x/N.y,dT=0)+on(b1,du2=0,dT=0);
u_ite1+=du1;
u_ite2+=du2;
p_ite+=dp;
T_ite+=dT;
k_ite+=dk;
ep_ite+=dep;
}
Upx=u_ite1;
Upy=u_ite2;
pp=p_ite;
Tp=T_ite;
kp=k_ite;
epp=ep_ite;
plot([Upx, Upy], p, coef=0.2, cmm=" [Ux, Uy] - p", WindowIndex=1);
``````

}

// Check
if (iter >= nbiter) break;
}
cout << "Total Time = " << clock()-T0 << endl;
how to do this?