Index: /issm/workshop/10_Mesh_generation/Square.exp
===================================================================
--- /issm/workshop/10_Mesh_generation/Square.exp	(revision 11036)
+++ /issm/workshop/10_Mesh_generation/Square.exp	(revision 11036)
@@ -0,0 +1,10 @@
+## Name:domainoutline
+## Icon:0
+# Points Count  Value
+5 1.
+# X pos Y pos
+0 0
+1 0
+1 1
+0 1
+0 0
Index: /issm/workshop/10_Mesh_generation/circles.m
===================================================================
--- /issm/workshop/10_Mesh_generation/circles.m	(revision 11036)
+++ /issm/workshop/10_Mesh_generation/circles.m	(revision 11036)
@@ -0,0 +1,7 @@
+function vel=circles(x,y),
+
+	u=4*x-2; v=4*y-2;
+
+	vel=tanh(30*(u.^2+v.^2-0.25)) ...
+		+tanh(30*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u-0.75).^2+(v+0.75).^2-0.25)) ...
+		+tanh(30*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u+0.75).^2+(v+0.75).^2-0.25)) ;
Index: /issm/workshop/12_Inverse_methods/Front.exp
===================================================================
--- /issm/workshop/12_Inverse_methods/Front.exp	(revision 11036)
+++ /issm/workshop/12_Inverse_methods/Front.exp	(revision 11036)
@@ -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/12_Inverse_methods/Square.par
===================================================================
--- /issm/workshop/12_Inverse_methods/Square.par	(revision 11036)
+++ /issm/workshop/12_Inverse_methods/Square.par	(revision 11036)
@@ -0,0 +1,31 @@
+%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.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness;
+md.geometry.surface=md.geometry.bed+md.geometry.thickness;
+
+disp('      creating drag');
+pos=find(md.mask.elementonfloatingice);
+md.friction.coefficient=200*ones(md.mesh.numberofvertices,1);
+md.friction.coefficient(md.mesh.elements(pos,:))=0;
+md.friction.p=ones(md.mesh.numberofelements,1);
+md.friction.q=ones(md.mesh.numberofelements,1);
+
+disp('      initial velocity');
+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.vel=zeros(md.mesh.numberofvertices,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/12_Inverse_methods/runme.m
===================================================================
--- /issm/workshop/12_Inverse_methods/runme.m	(revision 11036)
+++ /issm/workshop/12_Inverse_methods/runme.m	(revision 11036)
@@ -0,0 +1,51 @@
+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,'macayeal','all');
+	md.cluster=generic('name',oshostname,'np',2);
+	md=solve(md,DiagnosticSolutionEnum);
+	plotmodel(md,'data',md.materials.rheology_B,'caxis',[ 1.3 1.9]*10^8,'figure',1);
+	plotmodel(md,'data',md.results.DiagnosticSolution.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.vx_obs=md.results.DiagnosticSolution.Vx;
+	md.inversion.vy_obs=md.results.DiagnosticSolution.Vy;
+	md.inversion.vel_obs=md.results.DiagnosticSolution.Vel;
+
+	md=solve(md,DiagnosticSolutionEnum);
+	plotmodel(md,'data',md.materials.rheology_B,'caxis',[ 1.3 1.9]*10^8,'figure',1);
+	plotmodel(md,'data',md.results.DiagnosticSolution.Vel,'figure',2);  
+	save model2 md
+end
+if step==3
+	%invert for rheology
+	loadmodel('model2.mat');
+
+	%Set up inversion parameters
+	nsteps=40;
+	md.inversion.iscontrol=1;
+	md.inversion.control_parameters={'MaterialsRheologyBbar'};
+	md.inversion.nsteps=nsteps;
+	md.inversion.cost_functions=101*ones(nsteps,1);
+	md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,1);
+	md.inversion.maxiter_per_step=10*ones(nsteps,1);
+	md.inversion.step_threshold=.8*ones(nsteps,1);
+	md.inversion.gradient_scaling=10^7*ones(nsteps,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=solve(md,DiagnosticSolutionEnum);
+	plotmodel(md,'data',md.results.DiagnosticSolution.MaterialsRheologyBbar,'caxis',[ 1.3 1.9]*10^8,'figure',1);
+	plotmodel(md,'data',md.results.DiagnosticSolution.Vel,'figure',2);  
+end