gf_workspace('clear all');

L=100;H=100;
m=gfMesh('regular simplices',0:10:L, 0:10:H);

mfd=gfMeshFem(m,1);
mfd.set('fem',gfFem('FEM_PK(2,1)'));
mim=gf_mesh_im(m,gf_integ('IM_TRIANGLE(3)'));

pts=m.pts;
top_pid=find(abs(pts(2,:)-L)<=1e-6);
bot_pid=find(abs(pts(2,:)-0)<=1e-6);

top =gf_mesh_get(m,'faces from pid',top_pid);
bot =gf_mesh_get(m,'faces from pid',bot_pid);
all=gf_mesh_get(m,'outer faces');

m.set('region',101,top);
m.set('region',102,bot);
m.set('region',103,setdiff(all',union(top',bot','rows'),'rows')');

md=gf_model('real');
gf_model_set(md,'add fem variable','V',mfd);
gf_model_set(md,'add Laplacian brick',mim,'V');

gf_model_set(md,'add initialized data','V0',0);
gf_model_set(md,'add initialized data','V1',1);
gf_model_set(md,'add initialized data','Q0',0);

disp('1')
gf_model_set(md,'add Dirichlet condition with multipliers',mim,'V',mfd,101,'V0');
disp('2')
gf_model_set(md,'add Dirichlet condition with multipliers',mim,'V',mfd,102,'V1');
disp('3')
gf_model_set(md,'add normal source term brick',mim,'V','Q0',103);

disp('solve')
gf_model_get(md,'solve');
