Index: /issm/workshop/05Inversion/DomainOutline.exp
===================================================================
--- /issm/workshop/05Inversion/DomainOutline.exp	(revision 18051)
+++ /issm/workshop/05Inversion/DomainOutline.exp	(revision 18051)
@@ -0,0 +1,10 @@
+## Name:DomainOutline
+## Icon:0
+# Points Count  Value
+5 1.000000
+# X pos Y pos
+0 0
+1000000 0
+1000000 1000000
+0 1000000
+0 0
Index: /issm/workshop/05Inversion/Front.exp
===================================================================
--- /issm/workshop/05Inversion/Front.exp	(revision 18051)
+++ /issm/workshop/05Inversion/Front.exp	(revision 18051)
@@ -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/05Inversion/Square.par
===================================================================
--- /issm/workshop/05Inversion/Square.par	(revision 18051)
+++ /issm/workshop/05Inversion/Square.par	(revision 18051)
@@ -0,0 +1,24 @@
+%Start defining model parameters here
+
+disp('      creating thickness');
+hmin = 300;
+hmax = 1000;
+ymin = min(md.mesh.y);
+ymax = max(md.mesh.y);
+md.geometry.thickness = hmax+(hmin-hmax)*(md.mesh.y-ymin)/(ymax-ymin);
+md.geometry.base      = -md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness;
+md.geometry.surface   = md.geometry.base+md.geometry.thickness;
+
+disp('      creating drag');
+md.friction.coefficient=200*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 flow law paramter');
+md.materials.rheology_B=1.8*10^8*ones(md.mesh.numberofvertices,1);
+md.materials.rheology_B(find(md.mesh.x<md.mesh.y))=1.4*10^8;
+md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+
+disp('      creating boundary conditions');
+md=SetIceShelfBC(md,'Front.exp');
Index: /issm/workshop/05Inversion/runme.m
===================================================================
--- /issm/workshop/05Inversion/runme.m	(revision 18051)
+++ /issm/workshop/05Inversion/runme.m	(revision 18051)
@@ -0,0 +1,50 @@
+step=1;
+if step==1
+	%Generate observation
+	md = model;
+	md = triangle(md,'DomainOutline.exp',100000);
+	md = setmask(md,'all','');
+	md = parameterize(md,'Square.par');
+	md = setflowequation(md,'SSA','all');
+	md.cluster = generic('np',2);
+	md = solve(md,StressbalanceSolutionEnum);
+	plotmodel(md,'data',md.materials.rheology_B,'caxis',[ 1.3 1.9]*10^8,'figure',1);
+	plotmodel(md,'data',md.results.StressbalanceSolution.Vel,'figure',2);  
+	save model1 md
+end
+if step==2
+	%Modify rheology, now constant
+	loadmodel('model1.mat');
+	md.materials.rheology_B(:) = 1.8*10^8;
+
+	%results of previous run are taken as observations
+	md.inversion=m1qn3inversion();
+	md.inversion.vx_obs  = md.results.StressbalanceSolution.Vx;
+	md.inversion.vy_obs  = md.results.StressbalanceSolution.Vy;
+	md.inversion.vel_obs = md.results.StressbalanceSolution.Vel;
+
+	md = solve(md,StressbalanceSolutionEnum);
+	plotmodel(md,'data',md.materials.rheology_B,'caxis',[ 1.3 1.9]*10^8,'figure',1);
+	plotmodel(md,'data',md.results.StressbalanceSolution.Vel,'figure',2);  
+	save model2 md
+end
+if step==3
+	%invert for ice rigidity
+	loadmodel('model2.mat');
+
+	%Set up inversion parameters
+	maxsteps = 20;
+	md.inversion.iscontrol = 1;
+	md.inversion.control_parameters = {'MaterialsRheologyBbar'};
+	md.inversion.maxsteps = maxsteps;
+	md.inversion.cost_functions = [101];
+	md.inversion.cost_functions_coefficients = ones(md.mesh.numberofvertices,1);
+	md.inversion.min_parameters    = paterson(273)*ones(md.mesh.numberofvertices,1);
+	md.inversion.max_parameters    = paterson(200)*ones(md.mesh.numberofvertices,1);
+
+	%Go solve!
+	md.verbose=verbose(0);
+	md=solve(md,StressbalanceSolutionEnum);
+	plotmodel(md,'data',md.results.StressbalanceSolution.MaterialsRheologyBbar,'caxis',[ 1.3 1.9]*10^8,'figure',1);
+	plotmodel(md,'data',md.results.StressbalanceSolution.Vel,'figure',2);  
+end