Index: /issm/workshop/04IceflowModels/EISMINT.par
===================================================================
--- /issm/workshop/04IceflowModels/EISMINT.par	(revision 18047)
+++ /issm/workshop/04IceflowModels/EISMINT.par	(revision 18047)
@@ -0,0 +1,53 @@
+disp('      creating thickness');
+hmin=0.01;
+hmax=2756.7;
+radius=sqrt((md.mesh.x).^2+(md.mesh.y).^2);
+radiusmax=max(radius);
+radius(find(radius>(1.-10^-9)*radiusmax))=radiusmax;    %eliminate roundoff issues in next statement
+md.geometry.thickness=hmin*ones(size(md.mesh.x,1),1)+hmax*(4.*((1./2.)^(4./3.)*ones(size(md.mesh.x,1),1)-((radius)./(2.*radiusmax)).^(4./3.))).^(3./8.);
+md.geometry.base=0.*md.geometry.thickness;
+md.geometry.surface=md.geometry.base+md.geometry.thickness;
+
+disp('      creating drag');
+md.friction.coefficient=20.*ones(md.mesh.numberofvertices,1);
+md.friction.coefficient(find(md.mask.groundedice_levelset<0.))=0.;
+md.friction.p=ones(md.mesh.numberofelements,1);
+md.friction.q=ones(md.mesh.numberofelements,1);
+
+disp('      creating temperatures');
+tmin=238.15; %K
+st=1.67*10^-2/1000.; %k/m
+md.initialization.temperature=tmin+st*radius;
+md.basalforcings.geothermalflux=4.2*10^-2*ones(md.mesh.numberofvertices,1);
+
+disp('      creating flow law parameter');
+md.materials.rheology_B=6.81*10^7*ones(md.mesh.numberofvertices,1); %to have the same B as the analytical solution 
+md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1);
+
+disp('      creating surface mass balance');
+smb_max=0.5; %m/yr
+sb=10^-2/1000.; %m/yr/m
+rel=450.*1000.; %m
+md.surfaceforcings.mass_balance=min(smb_max,sb*(rel-radius));
+
+disp('      creating velocities');
+constant=0.3;
+md.inversion.vx_obs=constant/2.*md.mesh.x.*(md.geometry.thickness).^-1;
+md.inversion.vy_obs=constant/2.*md.mesh.y.*(md.geometry.thickness).^-1;
+md.inversion.vel_obs=sqrt((md.inversion.vx_obs).^2+(md.inversion.vy_obs).^2);
+md.initialization.vx=zeros(md.mesh.numberofvertices,1);
+md.initialization.vy=zeros(md.mesh.numberofvertices,1);
+md.initialization.vz=zeros(md.mesh.numberofvertices,1);
+md.initialization.pressure=zeros(md.mesh.numberofvertices,1);
+
+%Deal with boundary conditions:
+disp('      boundary conditions for stressbalance model:');
+md=SetMarineIceSheetBC(md,'./RoundFrontEISMINT.exp');
+
+radius=sqrt((md.mesh.x).^2+(md.mesh.y).^2);
+pos=find(radius==min(radius));
+md.mesh.x(pos)=0.; md.mesh.y(pos)=0.; %the closest node to the center is changed to be exactly at the center
+
+md.stressbalance.spcvx(pos)=0.;
+md.stressbalance.spcvy(pos)=0.;
+md.stressbalance.spcvz(pos)=0.;
Index: /issm/workshop/04IceflowModels/RoundFront.exp
===================================================================
--- /issm/workshop/04IceflowModels/RoundFront.exp	(revision 18047)
+++ /issm/workshop/04IceflowModels/RoundFront.exp	(revision 18047)
@@ -0,0 +1,47 @@
+## Name:
+## Icon:0
+# Points Count Value
+41 1
+# X pos Y pos
+1000000.000000 0.000000
+987688.340595 156434.465040
+951056.516295 309016.994375
+891006.524188 453990.499740
+809016.994375 587785.252292
+707106.781187 707106.781187
+587785.252292 809016.994375
+453990.499740 891006.524188
+309016.994375 951056.516295
+156434.465040 987688.340595
+0.000000 1000000.000000
+-156434.465040 987688.340595
+-309016.994375 951056.516295
+-453990.499740 891006.524188
+-587785.252292 809016.994375
+-707106.781187 707106.781187
+-809016.994375 587785.252292
+-891006.524188 453990.499740
+-951056.516295 309016.994375
+-987688.340595 156434.465040
+-1000000.000000 0.000000
+-987688.340595 -156434.465040
+-951056.516295 -309016.994375
+-891006.524188 -453990.499740
+-809016.994375 -587785.252292
+-707106.781187 -707106.781187
+-587785.252292 -809016.994375
+-453990.499740 -891006.524188
+-309016.994375 -951056.516295
+-156434.465040 -987688.340595
+-0.000000 -1000000.000000
+156434.465040 -987688.340595
+309016.994375 -951056.516295
+453990.499740 -891006.524188
+587785.252292 -809016.994375
+707106.781187 -707106.781187
+809016.994375 -587785.252292
+891006.524188 -453990.499740
+951056.516295 -309016.994375
+987688.340595 -156434.465040
+1000000.000000 0.000000
+
Index: /issm/workshop/04IceflowModels/Square.exp
===================================================================
--- /issm/workshop/04IceflowModels/Square.exp	(revision 18047)
+++ /issm/workshop/04IceflowModels/Square.exp	(revision 18047)
@@ -0,0 +1,10 @@
+## Name:domainoutline
+## Icon:0
+# Points Count  Value
+5 1.
+# X pos Y pos
+0 0
+1000000 0
+1000000 1000000
+0 1000000
+0 0
Index: /issm/workshop/04IceflowModels/SquareFront.exp
===================================================================
--- /issm/workshop/04IceflowModels/SquareFront.exp	(revision 18047)
+++ /issm/workshop/04IceflowModels/SquareFront.exp	(revision 18047)
@@ -0,0 +1,10 @@
+## Name:icefront
+## Icon:0
+# Points Count  Value
+5 1.
+# X pos Y pos
+-1000 900000
+-1000 1100000
+1100000 1100000
+1100000 900000
+-1000 900000
Index: /issm/workshop/04IceflowModels/SquareShelf.par
===================================================================
--- /issm/workshop/04IceflowModels/SquareShelf.par	(revision 18047)
+++ /issm/workshop/04IceflowModels/SquareShelf.par	(revision 18047)
@@ -0,0 +1,55 @@
+%Start defining model parameters here
+
+%Geometry
+hmin=300;
+hmax=1000;
+ymin=min(md.mesh.y);
+ymax=max(md.mesh.y);
+xmin=min(md.mesh.x);
+xmax=max(md.mesh.x);
+md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y-ymin)/(ymax-ymin)+0.1*(hmin-hmax)*(md.mesh.x-xmin)/(xmax-xmin);
+md.geometry.base=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness;
+md.geometry.surface=md.geometry.base+md.geometry.thickness;
+
+%Initial velocity and pressure
+x     = transpose(ncread('../Data/SquareShelf.nc','x'));
+y     = transpose(ncread('../Data/SquareShelf.nc','y'));
+vx    = transpose(ncread('../Data/SquareShelf.nc','vx'));
+vy    = transpose(ncread('../Data/SquareShelf.nc','vy'));
+index = transpose(ncread('../Data/SquareShelf.nc','index'));
+md.initialization.vx=InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y);
+md.initialization.vy=InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y);
+clear vx vy x y index;
+md.initialization.vz=zeros(md.mesh.numberofvertices,1);
+md.initialization.pressure=zeros(md.mesh.numberofvertices,1);
+
+%Materials
+md.initialization.temperature=(273-20)*ones(md.mesh.numberofvertices,1);
+md.materials.rheology_B=paterson(md.initialization.temperature);
+md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+
+%Friction
+md.friction.coefficient=20*ones(md.mesh.numberofvertices,1);
+md.friction.coefficient(find(md.mask.groundedice_levelset<0.))=0.;
+md.friction.p=ones(md.mesh.numberofelements,1);
+md.friction.q=ones(md.mesh.numberofelements,1);
+
+%Numerical parameters
+md.stressbalance.viscosity_overshoot=0.3;
+md.masstransport.stabilization=1;
+md.thermal.stabilization=1;
+md.settings.waitonlock=30;
+md.verbose=verbose(0);
+md.stressbalance.restol=0.10;
+md.steadystate.reltol=0.02;
+md.stressbalance.reltol=0.02;
+md.stressbalance.abstol=NaN;
+md.timestepping.time_step=1;
+md.timestepping.final_time=3;
+
+%Boundary conditions:
+md=SetIceShelfBC(md,'../Exp/SquareFront.exp');
+
+%Change name so that no test have the same name
+A=dbstack;
+if (length(A)>2), md.miscellaneous.name=A(3).file(1:end-2); end
Index: /issm/workshop/04IceflowModels/eismint.m
===================================================================
--- /issm/workshop/04IceflowModels/eismint.m	(revision 18047)
+++ /issm/workshop/04IceflowModels/eismint.m	(revision 18047)
@@ -0,0 +1,51 @@
+md=model();
+
+%Create mesh with roundmesh
+md=roundmesh(md,750000,30000);
+
+%Set mask
+md=setmask(md,'','');
+
+%Parameterize model
+md=parameterize(md,'EISMINT.par');
+
+%We extrude the model to have a 3d model
+md=extrude(md,10,1);
+
+%Set ice flow approximation
+...
+
+%Create boundary conditions: zero velocity on the bed
+pos=find(md.mesh.vertexonbed);
+md.stressbalance.spcvx(pos)=0;
+md.stressbalance.spcvy(pos)=0;
+md.stressbalance.spcvz(pos)=0;
+
+%Go Solve
+md.cluster=generic('np',2);
+md.verbose.convergence=1;
+md=solve(md,StressbalanceSolutionEnum());
+vel=DepthAverage(md,sqrt(md.results.StressbalanceSolution.Vx.^2+md.results.StressbalanceSolution.Vy.^2));
+
+%Calculate analytical velocity
+constant=0.3;
+vx_obs=constant/2*md.mesh.x.*(md.geometry.thickness).^-1;
+vy_obs=constant/2*md.mesh.y.*(md.geometry.thickness).^-1;
+vel_obs=sqrt(vx_obs.^2+vy_obs.^2);
+vel_obs=project2d(md,vel_obs,1);
+
+plotmodel(md,...
+	'data',vel    ,'view',2,'caxis',[0 200],'title','Modelled velocity',...
+	'data',vel_obs,'view',2,'caxis',[0 200],'title','Analytical velocity',...
+	'data',abs(vel-vel_obs)./(vel_obs+eps)*100,'caxis',[0 30],'view',2,'title','Relative misfit (%)');
+
+subplot(2,2,4)
+hold on;
+plot(sqrt((md.mesh.x2d).^2+(md.mesh.y2d).^2),vel,'r.');
+plot(sqrt((md.mesh.x2d).^2+(md.mesh.y2d).^2),vel_obs,'b.');
+title('Analytical vs calculated velocity');
+xlabel('distance to the center of the icesheet [m]');
+ylabel('velocity (m/yr)');
+legend('calculated velocity','exact velocity');
+axis([0 750000 0 200]);
+hold off;
Index: /issm/workshop/04IceflowModels/runme.m
===================================================================
--- /issm/workshop/04IceflowModels/runme.m	(revision 18047)
+++ /issm/workshop/04IceflowModels/runme.m	(revision 18047)
@@ -0,0 +1,7 @@
+md=triangle(model(),'Square.exp',50000.);
+md=setmask(md,'all','');
+md=parameterize(md,'SquareShelf.par');
+md=extrude(md,3,1);
+%Set flow equation
+...
+%Solve