Index: /issm/trunk/src/m/utils/Qmu/lrel_mmf.m
===================================================================
--- /issm/trunk/src/m/utils/Qmu/lrel_mmf.m	(revision 4701)
+++ /issm/trunk/src/m/utils/Qmu/lrel_mmf.m	(revision 4701)
@@ -0,0 +1,96 @@
+%  set up a local reliability study, like might be done in Pig.par
+
+%%  a variety of variables
+
+%  seems to be a Matlab bug here (on Linux, not WinXP) -- unless
+%  the class has been called, "empty" method can not be found
+normal_uncertain;
+continuous_design;
+continuous_state;
+linear_inequality_constraint;
+linear_equality_constraint;
+response_function;
+objective_function;
+least_squares_term;
+nonlinear_inequality_constraint;
+nonlinear_equality_constraint;
+
+md.variables=struct();
+md.variables.nuv=normal_uncertain.empty();
+%md.variables.nuv(end+1)=normal_uncertain('rho_ice',917,45.85);
+%md.variables.nuv(end+1)=normal_uncertain('rho_water',1023,51.15);
+%md.variables.nuv(end+1)=normal_uncertain('heatcapacity',2009,100.45);
+%md.variables.nuv(end+1)=normal_uncertain('thermalconductivity',2.2,0.11);
+%md.variables.nuv(end+1)=normal_uncertain('gravity',9.8,0.49);
+md.variables.nuv(end+1)=normal_uncertain('thickness',1,0.05);
+%md.variables.nuv(end+1)=normal_uncertain('drag',1,0.05);
+
+%%  a variety of responses
+
+md.responses=struct();
+md.responses.rf =response_function.empty();
+md.responses.rf (end+1)=response_function('min_vx',[],[0.0001 0.001 0.01 0.25 0.5 0.75 0.99 0.999 0.9999]);
+md.responses.rf (end+1)=response_function('max_vx',[],[0.0001 0.001 0.01 0.25 0.5 0.75 0.99 0.999 0.9999]);
+md.responses.rf (end+1)=response_function('max_abs_vx',[],[0.0001 0.001 0.01 0.25 0.5 0.75 0.99 0.999 0.9999]);
+md.responses.rf (end+1)=response_function('min_vy',[],[0.0001 0.001 0.01 0.25 0.5 0.75 0.99 0.999 0.9999]);
+md.responses.rf (end+1)=response_function('max_vy',[],[0.0001 0.001 0.01 0.25 0.5 0.75 0.99 0.999 0.9999]);
+md.responses.rf (end+1)=response_function('max_abs_vy',[],[0.0001 0.001 0.01 0.25 0.5 0.75 0.99 0.999 0.9999]);
+md.responses.rf (end+1)=response_function('min_vel',[],[0.0001 0.001 0.01 0.25 0.5 0.75 0.99 0.999 0.9999]);
+md.responses.rf (end+1)=response_function('max_vel',[],[0.0001 0.001 0.01 0.25 0.5 0.75 0.99 0.999 0.9999]);
+md.responses.rf (end+1)=response_function('mass_flux1',[],[0.0001 0.001 0.01 0.25 0.5 0.75 0.99 0.999 0.9999]);
+md.responses.rf (end+1)=response_function('mass_flux_2',[],[0.0001 0.001 0.01 0.25 0.5 0.75 0.99 0.999 0.9999]);
+md.responses.rf (end+1)=response_function('mass_flux(3)',[],[0.0001 0.001 0.01 0.25 0.5 0.75 0.99 0.999 0.9999]);
+md.responses.rf (end+1)=response_function('mass_flux4 (repeat)',[],[0.0001 0.001 0.01 0.25 0.5 0.75 0.99 0.999 0.9999]);
+md.responses.rf (end+1)=response_function('mass_flux-5',[],[0.0001 0.001 0.01 0.25 0.5 0.75 0.99 0.999 0.9999]);
+md.responses.rf (end+1)=response_function('mass_flux^6',[],[0.0001 0.001 0.01 0.25 0.5 0.75 0.99 0.999 0.9999]);
+md.responses.rf (end+1)=response_function('mass_flux[7]',[],[0.0001 0.001 0.01 0.25 0.5 0.75 0.99 0.999 0.9999]);
+
+%%  create mass flux profile
+
+%plotmodel(md,'data','mesh')
+%expcreateprofile('mass_flux.exp')
+%expdisp('mass_flux.exp')
+%md.qmu_mass_flux_profile='mass_flux.exp';
+%md.qmu_mass_flux_profile={'mass_flux.exp','mass_flux2.exp','mass_flux3.exp'};
+md.qmu_mass_flux_profile={'mass_flux.exp','mass_flux2.exp','mass_flux3.exp','mass_flux.exp','mass_flux5.exp','mass_flux6.exp','mass_flux7.exp'};
+
+%%  nond_local_reliability study
+
+md.qmu_method     =dakota_method('nond_l');
+
+%%  a variety of parameters
+
+%md.qmu_params.evaluation_concurrency=4;
+md.qmu_params.evaluation_concurrency=1;
+md.qmu_params.analysis_driver='';
+md.qmu_params.analysis_components='';
+md.qmu_params.interval_type='forward';
+md.qmu_params.fd_gradient_step_size=0.01;
+
+md.qmu_analysis=1;
+md.npart=10;
+if isempty(md.adjacency)
+	md=adjacency(md);
+end
+if isempty(md.part)
+%    md.part=partitioner(md,'package','metis','npart',md.npart);
+    md.part=partitioner(md,'package','chaco','npart',md.npart,'weighting','on');
+% SpawnCore.m assumes partition vector starting at zero
+    md.part=md.part-1;
+end
+md.eps_rel=1.e-5;
+md.verbose=1;
+md.cluster='none';
+
+md.qmu
+
+%%  sample analysis
+
+%md=solve(md,'analysis_type','diagnostic');
+
+%plotmodel(md,'data','mesh')
+%plotmodel(md,'data',md.part)
+%plotmodel(md,'data','mesh','partitionedges','on','linewidth',2)
+%part_hist(md.part,md.vwgt)
+%plotmodel(md,'data',log10(md.results.dakota.dresp_out(9).impfac(md.part+1)))
+