b:=10;      /* total length of YSZ                       */
h:=4;       /* thickness of YSZ                          */
w:=2;       /* width of working electrode               */
a:=4;       /* particle radius                           */
vw:=1;      /* potential at working electrode             */
vc:=0;      /* potential at counter electrode             */
null:=0.000001;     /* null condition                       */
aa:=20;        /* gridpoints per boundary                 */

/*    |-----------------b-----------------|
                                 
                     a,h
                      | (a+w),h
      |------a--------|   |
     null,h           |-w-|              b,h
---   |---------------*****---------------|    potential: vw
 |    |   (2)          (3)       WE  (1)  |
 h    |                                   |
 |    |(2)                             (1)|
 |    |                (4)                |
---   |***********************************|    potential: vc
     null,0           CE                 b,0
*/
/* rectangle between (0,0) and (b,h) */

border(4, null, b,30)           begin x:=t; y:=0; end;
border(1, 0, h,   30)           begin x:=b; y:=t; end;
border(1, b, a+w, 30)           begin x:=t; y:=h; end;
border(3, a+w,a,  30)           begin x:=t; y:=h; end;
border(2, a,null, 30)           begin x:=t; y:=h; end;
border(2, h, 0,   30)           begin x:=null; y:=t; end;

buildmesh(2000);

f=x;
g=y;
savemesh('constr1.mes');
invx=1/x;
solve(u)
begin
  onbdy(1) dnu(u)=0;
  onbdy(2) dnu(u)=0;
  onbdy(3) u=vw;
  onbdy(4) u=vc;
/*  pde(u) laplace(u) = 0.;*/
  pde(u) laplace(u) + dx(u)*1/x = 0.;
end; 
/*adaptmesh(u);*/

plot(u);

/* solve adjoint problem to display current lines */

jx = dx(u);
jy = dy(u);

solve(v)
begin
  onbdy(1) v = 0;
  onbdy(2) v = 1;
  onbdy(3) dnu(v)=0;
  onbdy(4) dnu(v)=0;
  pde(v) laplace(v) = 0.;
end;

plot(v);
save('constr1.out',u);