Index: /issm/trunk/src/m/classes/@continuous_design/continuous_design.m =================================================================== --- /issm/trunk/src/m/classes/@continuous_design/continuous_design.m (revision 26) +++ /issm/trunk/src/m/classes/@continuous_design/continuous_design.m (revision 27) @@ -7,7 +7,9 @@ properties descriptor=''; - initpt = NaN; - lower = NaN; - upper = NaN; + initpt = 0.; + lower =-Inf; + upper = Inf; + scale_type='none'; + scale = 1.; end @@ -35,13 +37,22 @@ otherwise cdv.descriptor=varargin{1}; - cdv.initpt =varargin{2}; - if (nargin >= 3) - cdv.lower =varargin{3}; - if (nargin >= 4) - cdv.upper =varargin{4}; - if (nargin > 4) - warning('continuous_design:extra_arg',... - 'Extra arguments for object of class ''%s''.',... - class(cdv)); + + if (nargin >= 2) + cdv.initpt =varargin{2}; + if (nargin >= 3) + cdv.lower =varargin{3}; + if (nargin >= 4) + cdv.upper =varargin{4}; + if (nargin >= 5) + cdv.scale_type=varargin{5}; + if (nargin >= 6) + cdv.scale =varargin{6}; + if (nargin > 6) + warning('continuous_design:extra_arg',... + 'Extra arguments for object of class ''%s''.',... + class(cdv)); + end + end + end end end @@ -50,5 +61,5 @@ end - function [desc] =dvar_desc(cdv) + function [desc] =prop_desc(cdv) desc=cell(size(cdv)); for i=1:numel(cdv) @@ -57,32 +68,47 @@ desc=allempty(desc); end - function [initpt]=dvar_initpt(cdv) + function [initpt]=prop_initpt(cdv) initpt=zeros(size(cdv)); for i=1:numel(cdv) initpt(i)=cdv(i).initpt; end + initpt=allequal(initpt,0.); end - function [lower] =dvar_lower(cdv) + function [lower] =prop_lower(cdv) lower=zeros(size(cdv)); for i=1:numel(cdv) lower(i)=cdv(i).lower; end - lower=allnan(lower); + lower=allequal(lower,-Inf); end - function [upper] =dvar_upper(cdv) + function [upper] =prop_upper(cdv) upper=zeros(size(cdv)); for i=1:numel(cdv) upper(i)=cdv(i).upper; end - upper=allnan(upper); + upper=allequal(upper, Inf); end - function [mean] =dvar_mean(cdv) + function [mean] =prop_mean(cdv) mean=[]; end - function [stddev]=dvar_stddev(cdv) + function [stddev]=prop_stddev(cdv) stddev=[]; end - function [initst]=dvar_initst(cdv) + function [initst]=prop_initst(cdv) initst=[]; + end + function [stype] =prop_stype(cdv) + stype=cell(size(cdv)); + for i=1:numel(cdv) + stype(i)=cellstr(cdv(i).scale_type); + end + stype=allequal(stype,'none'); + end + function [scale] =prop_scale(cdv) + scale=zeros(size(cdv)); + for i=1:numel(cdv) + scale(i)=cdv(i).scale; + end + scale=allequal(scale,1.); end end Index: /issm/trunk/src/m/classes/@continuous_design/display.m =================================================================== --- /issm/trunk/src/m/classes/@continuous_design/display.m (revision 26) +++ /issm/trunk/src/m/classes/@continuous_design/display.m (revision 27) @@ -20,5 +20,7 @@ disp(sprintf(' initpt: %g' ,cdv(i).initpt)); disp(sprintf(' lower: %g' ,cdv(i).lower)); - disp(sprintf(' upper: %g\n' ,cdv(i).upper)); + disp(sprintf(' upper: %g' ,cdv(i).upper)); + disp(sprintf(' scale_type: ''%s''' ,cdv(i).scale_type)); + disp(sprintf(' scale: %g' ,cdv(i).scale)); end Index: /issm/trunk/src/m/classes/@continuous_state/continuous_state.m =================================================================== --- /issm/trunk/src/m/classes/@continuous_state/continuous_state.m (revision 26) +++ /issm/trunk/src/m/classes/@continuous_state/continuous_state.m (revision 27) @@ -7,7 +7,7 @@ properties descriptor=''; - initst = NaN; - lower = NaN; - upper = NaN; + initst = 0.; + lower =-Inf; + upper = Inf; end @@ -35,13 +35,16 @@ otherwise csv.descriptor=varargin{1}; - csv.initst =varargin{2}; - if (nargin >= 3) - csv.lower =varargin{3}; - if (nargin >= 4) - csv.upper =varargin{4}; - if (nargin > 4) - warning('continuous_state:extra_arg',... - 'Extra arguments for object of class ''%s''.',... - class(csv)); + + if (nargin >= 2) + csv.initst =varargin{2}; + if (nargin >= 3) + csv.lower =varargin{3}; + if (nargin >= 4) + csv.upper =varargin{4}; + if (nargin > 4) + warning('continuous_state:extra_arg',... + 'Extra arguments for object of class ''%s''.',... + class(csv)); + end end end @@ -50,5 +53,5 @@ end - function [desc] =dvar_desc(csv) + function [desc] =prop_desc(csv) desc=cell(size(csv)); for i=1:numel(csv) @@ -57,32 +60,39 @@ desc=allempty(desc); end - function [initpt]=dvar_initpt(csv) + function [initpt]=prop_initpt(csv) initpt=[]; end - function [lower] =dvar_lower(csv) + function [lower] =prop_lower(csv) lower=zeros(size(csv)); for i=1:numel(csv) lower(i)=csv(i).lower; end - lower=allnan(lower); + lower=allequal(lower,-Inf); end - function [upper] =dvar_upper(csv) + function [upper] =prop_upper(csv) upper=zeros(size(csv)); for i=1:numel(csv) upper(i)=csv(i).upper; end - upper=allnan(upper); + upper=allequal(upper, Inf); end - function [mean] =dvar_mean(csv) + function [mean] =prop_mean(csv) mean=[]; end - function [stddev]=dvar_stddev(csv) + function [stddev]=prop_stddev(csv) stddev=[]; end - function [initst]=dvar_initst(csv) + function [initst]=prop_initst(csv) initst=zeros(size(csv)); for i=1:numel(csv) initst(i)=csv(i).initst; end + initst=allequal(initst,0.); + end + function [stype] =prop_stype(csv) + stype={}; + end + function [scale] =prop_scale(csv) + scale=[]; end end Index: /issm/trunk/src/m/classes/@dakota_method/dakota_method.m =================================================================== --- /issm/trunk/src/m/classes/@dakota_method/dakota_method.m (revision 27) +++ /issm/trunk/src/m/classes/@dakota_method/dakota_method.m (revision 27) @@ -0,0 +1,699 @@ +% +% definition for the dakota_method class. +% +% [dm]=dakota_method(method) +% +classdef dakota_method + properties + method =''; + type =''; + variables={}; + lcspec ={}; + responses={}; + ghspec ={}; + params =struct(); + end + + methods + function [dm]=dakota_method(method) + + switch nargin + +% create a default object + + case 0 + +% copy the object or create the object from the input + + case 1 + if (nargin == 1) && isa(method,'dakota_method') + dm=method; + else + mlist={... + 'dot_bfgs',... + 'dot_frcg',... + 'dot_mmfd',... + 'dot_slp',... + 'dot_sqp',... + 'npsol_sqp',... + 'conmin_frcg',... + 'conmin_mfd',... + 'optpp_cg',... + 'optpp_q_newton',... + 'optpp_fd_newton',... + 'optpp_newton',... + 'optpp_pds',... + 'asynch_pattern_search',... + 'coliny_cobyla',... + 'coliny_direct',... + 'coliny_ea',... + 'coliny_pattern_search',... + 'coliny_solis_wets',... + 'ncsu_direct',... + 'surrogate_based_local',... + 'surrogate_based_global',... + 'moga',... + 'soga',... + 'nl2sol',... + 'nlssol_sqp',... + 'optpp_g_newton',... + 'nond_sampling',... + 'nond_local_reliability',... + 'nond_global_reliability',... + 'nond_polynomial_chaos',... + 'nond_stoch_collocation',... + 'nond_evidence',... + 'dace',... + 'fsu_quasi_mc',... + 'fsu_cvt',... + 'vector_parameter_study',... + 'list_parameter_study',... + 'centered_parameter_study',... + 'multidim_parameter_study',... + }; + + mlist2={}; + for i=1:length(mlist) + if strncmpi(method,mlist{i},length(method)) + mlist2(end+1)=mlist(i); + end + end + +% check for a unique match in the list of methods + + switch length(mlist2) + case 0 + error('Unrecognized method: ''%s''.',... + method); + case 1 + dm.method=mlist2{1}; + otherwise + error('Non-unique method: ''%s'' matches %s.',... + method,string_cell(mlist2)); + end + +% assign the default values for the method + + switch dm.method + case {'dot_bfgs',... + 'dot_frcg'} + dm.type ='dot'; + dm.variables={'cdv','csv'}; + dm.lcspec ={}; + dm.responses={'of'}; + dm.ghspec ={'grad'}; + dm.params.max_iterations=false; + dm.params.max_function_evaluations=false; + dm.params.convergence_tolerance=false; + dm.params.constraint_tolerance=false; + dm.params.output=false; + dm.params.speculative=false; + dm.params.scaling=false; + dm.params.optimization_type='minimize'; + case {'dot_mmfd',... + 'dot_slp',... + 'dot_sqp'} + dm.type ='dot'; + dm.variables={'cdv','csv'}; + dm.lcspec ={'lic','lec'}; + dm.responses={'of','nic','nec'}; + dm.ghspec ={'grad'}; + dm.params.max_iterations=false; + dm.params.max_function_evaluations=false; + dm.params.convergence_tolerance=false; + dm.params.constraint_tolerance=false; + dm.params.output=false; + dm.params.speculative=false; + dm.params.scaling=false; + dm.params.optimization_type='minimize'; + + case {'npsol_sqp'} + dm.type ='npsol'; + dm.variables={'cdv','csv'}; + dm.lcspec ={'lic','lec'}; + dm.responses={'of','nic','nec'}; + dm.ghspec ={'grad'}; + dm.params.max_iterations=false; + dm.params.max_function_evaluations=false; + dm.params.convergence_tolerance=false; + dm.params.constraint_tolerance=false; + dm.params.output=false; + dm.params.speculative=false; + dm.params.scaling=false; + dm.params.verify_level=-1; + dm.params.function_precision=1.e-10; + dm.params.linesearch_tolerance=0.9; + + case {'conmin_frcg'} + dm.type ='conmin'; + dm.variables={'cdv','csv'}; + dm.lcspec ={}; + dm.responses={'of'}; + dm.ghspec ={'grad'}; + dm.params.max_iterations=false; + dm.params.max_function_evaluations=false; + dm.params.convergence_tolerance=false; + dm.params.constraint_tolerance=false; + dm.params.output=false; + dm.params.speculative=false; + dm.params.scaling=false; + case {'conmin_mfd'} + dm.type ='conmin'; + dm.variables={'cdv','csv'}; + dm.lcspec ={'lic','lec'}; + dm.responses={'of','nic','nec'}; + dm.ghspec ={'grad'}; + dm.params.max_iterations=false; + dm.params.max_function_evaluations=false; + dm.params.convergence_tolerance=false; + dm.params.constraint_tolerance=false; + dm.params.output=false; + dm.params.speculative=false; + dm.params.scaling=false; + + case {'optpp_cg'} + dm.type ='optpp'; + dm.variables={'cdv','csv'}; + dm.lcspec ={}; + dm.responses={'of'}; + dm.ghspec ={'grad'}; + dm.params.max_iterations=false; + dm.params.max_function_evaluations=false; + dm.params.convergence_tolerance=false; + dm.params.output=false; + dm.params.speculative=false; + dm.params.scaling=false; + dm.params.max_step=1000.; + dm.params.gradient_tolerance=1.e-4; + case {'optpp_q_newton',... + 'optpp_fd_newton',... + 'optpp_newton'} + dm.type ='optpp'; + dm.variables={'cdv','csv'}; + dm.lcspec ={'lic','lec'}; + dm.responses={'of','nic','nec'}; + dm.ghspec ={'grad'}; + dm.params.max_iterations=false; + dm.params.max_function_evaluations=false; + dm.params.convergence_tolerance=false; + dm.params.output=false; + dm.params.speculative=false; + dm.params.scaling=false; + dm.params.value_based_line_search=false; + dm.params.gradient_based_line_search=false; + dm.params.trust_region=false; + dm.params.tr_pds=false; + dm.params.max_step=1000.; + dm.params.gradient_tolerance=1.e-4; + dm.params.merit_function='argaez_tapia'; + dm.params.central_path=dm.params.merit_function; + dm.params.steplength_to_boundary=0.99995; + dm.params.centering_parameter=0.2; + case {'optpp_pds'} + dm.type ='optpp'; + dm.variables={'cdv','csv'}; + dm.lcspec ={}; + dm.responses={'of'}; + dm.ghspec ={'grad'}; + dm.params.max_iterations=false; + dm.params.max_function_evaluations=false; + dm.params.convergence_tolerance=false; + dm.params.output=false; + dm.params.speculative=false; + dm.params.scaling=false; + dm.params.search_scheme_size=32; + + case {'asynch_pattern_search'} + dm.type ='apps'; + dm.variables={'cdv','csv'}; + dm.lcspec ={'lic','lec'}; + dm.responses={'of','nic','nec'}; + dm.ghspec ={'grad'}; + dm.params.max_function_evaluations=false; + dm.params.constraint_tolerance=false; + dm.params.output=false; + dm.params.scaling=false; + dm.params.initial_delta=1.0; + dm.params.threshold_delta=0.01; + dm.params.contraction_factor=0.5; + dm.params.solution_target=false; + dm.params.synchronization='nonblocking'; + dm.params.merit_function='merit2_smooth'; + dm.params.constraint_penalty=1.0; + dm.params.smoothing_factor=1.0; + + case {'coliny_cobyla'} + dm.type ='coliny'; + dm.variables={'cdv','csv'}; + dm.lcspec ={}; + dm.responses={'of','nic','nec'}; + dm.ghspec ={'grad'}; + dm.params.max_iterations=false; + dm.params.max_function_evaluations=false; + dm.params.convergence_tolerance=false; + dm.params.output=false; + dm.params.scaling=false; + dm.params.show_misc_options=false; + dm.params.misc_options={}; + dm.params.solution_accuracy=-Inf; + dm.params.initial_delta=[]; + dm.params.threshold_delta=[]; + case {'coliny_direct'} + dm.type ='coliny'; + dm.variables={'cdv','csv'}; + dm.lcspec ={}; + dm.responses={'of','nic','nec'}; + dm.ghspec ={'grad'}; + dm.params.max_iterations=false; + dm.params.max_function_evaluations=false; + dm.params.convergence_tolerance=false; + dm.params.output=false; + dm.params.scaling=false; + dm.params.show_misc_options=false; + dm.params.misc_options={}; + dm.params.solution_accuracy=-Inf; + dm.params.division='major_dimension'; + dm.params.global_balance_parameter=0.0; + dm.params.local_balance_parameter=1.e-8; + dm.params.max_boxsize_limit=0.0; + dm.params.min_boxsize_limit=0.0001; + dm.params.constraint_penalty=1000.0; + case {'coliny_ea'} + dm.type ='coliny'; + dm.variables={'cdv','csv'}; + dm.lcspec ={}; + dm.responses={'of','nic','nec'}; + dm.ghspec ={'grad'}; + dm.params.max_iterations=false; + dm.params.max_function_evaluations=false; + dm.params.convergence_tolerance=false; + dm.params.output=false; + dm.params.scaling=false; + dm.params.show_misc_options=false; + dm.params.misc_options={}; + dm.params.solution_accuracy=-Inf; + dm.params.seed=false; + dm.params.population_size=50; + dm.params.initialization_type='unique_random'; + dm.params.fitness_type='linear_rank'; + dm.params.replacement_type='elitist'; + dm.params.random=[]; + dm.params.chc=[]; + dm.params.elitist=[]; + dm.params.new_solutions_generated='population_size - replacement_size'; + dm.params.crossover_type='two_point'; + dm.params.crossover_rate=0.8; + dm.params.mutation_type='offset_normal'; + dm.params.mutation_scale=0.1; + dm.params.mutation_range=1; + dm.params.dimension_ratio=1.0; + dm.params.mutation_rate=1.0; + dm.params.non_adaptive=false; + case {'coliny_pattern_search'} + dm.type ='coliny'; + dm.variables={'cdv','csv'}; + dm.lcspec ={}; + dm.responses={'of','nic','nec'}; + dm.ghspec ={'grad'}; + dm.params.max_iterations=false; + dm.params.max_function_evaluations=false; + dm.params.convergence_tolerance=false; + dm.params.output=false; + dm.params.scaling=false; + dm.params.show_misc_options=false; + dm.params.misc_options={}; + dm.params.solution_accuracy=-Inf; + dm.params.stochastic=false; + dm.params.seed=false; + dm.params.initial_delta=[]; + dm.params.threshold_delta=[]; + dm.params.constraint_penalty=1.0; + dm.params.constant_penalty=false; + dm.params.pattern_basis='coordinate'; + dm.params.total_pattern_size=false; + dm.params.no_expansion=false; + dm.params.expand_after_success=1; + dm.params.contraction_factor=0.5; + dm.params.synchronization='nonblocking'; + dm.params.exploratory_moves='basic_pattern'; + case {'coliny_solis_wets'} + dm.type ='coliny'; + dm.variables={'cdv','csv'}; + dm.lcspec ={}; + dm.responses={'of','nic','nec'}; + dm.ghspec ={'grad'}; + dm.params.max_iterations=false; + dm.params.max_function_evaluations=false; + dm.params.convergence_tolerance=false; + dm.params.output=false; + dm.params.scaling=false; + dm.params.show_misc_options=false; + dm.params.misc_options={}; + dm.params.solution_accuracy=-Inf; + dm.params.seed=false; + dm.params.initial_delta=[]; + dm.params.threshold_delta=[]; + dm.params.no_expansion=false; + dm.params.expand_after_success=5; + dm.params.contract_after_failure=3; + dm.params.contraction_factor=0.5; + dm.params.constraint_penalty=1.0; + dm.params.constant_penalty=false; + + case {'ncsu_direct'} + dm.type ='ncsu'; + dm.variables={'cdv','csv'}; + dm.lcspec ={'lic','lec'}; % ? + dm.responses={'of','nic','nec'}; % ? + dm.ghspec ={'grad'}; + dm.params.max_iterations=false; + dm.params.max_function_evaluations=false; + dm.params.scaling=false; + dm.params.solution_accuracy=0.; + dm.params.min_boxsize_limit=1.e-8; + dm.params.vol_boxsize_limit=1.e-8; + +% case {'surrogate_based_local',... +% 'surrogate_based_global'} + + case {'moga'} + dm.type ='jega'; + dm.variables={'cdv','csv'}; + dm.lcspec ={'lic','lec'}; + dm.responses={'of','nic','nec'}; + dm.ghspec ={'grad'}; + dm.params.max_iterations=false; + dm.params.max_function_evaluations=false; + dm.params.output=false; + dm.params.scaling=false; + dm.params.seed=false; + dm.params.log_file='JEGAGlobal.log'; + dm.params.population_size=50; + dm.params.print_each_pop=false; + dm.params.output='normal'; + dm.params.initialization_type='unique_random'; + dm.params.mutation_type='replace_uniform'; + dm.params.mutation_scale=0.15; + dm.params.mutation_rate=0.08; + dm.params.replacement_type=''; + dm.params.below_limit=6; + dm.params.shrinkage_percentage=0.9; + dm.params.crossover_type='shuffle_random'; + dm.params.multi_point_binary=[]; + dm.params.multi_point_parameterized_binary=[]; + dm.params.multi_point_real=[]; + dm.params.shuffle_random=[]; + dm.params.num_parents=2; + dm.params.num_offspring=2; + dm.params.crossover_rate=0.8; + dm.params.fitness_type=''; + dm.params.niching_type=false; + dm.params.radial=[0.01]; + dm.params.distance=[0.01]; + dm.params.metric_tracker=false; + dm.params.percent_change=0.1; + dm.params.num_generations=10; + dm.params.postprocessor_type=false; + dm.params.orthogonal_distance=[0.01]; + case {'soga'} + dm.type ='jega'; + dm.variables={'cdv','csv'}; + dm.lcspec ={'lic','lec'}; + dm.responses={'of','nic','nec'}; + dm.ghspec ={'grad'}; + dm.params.max_iterations=false; + dm.params.max_function_evaluations=false; + dm.params.output=false; + dm.params.scaling=false; + dm.params.seed=false; + dm.params.log_file='JEGAGlobal.log'; + dm.params.population_size=50; + dm.params.print_each_pop=false; + dm.params.output='normal'; + dm.params.initialization_type='unique_random'; + dm.params.mutation_type='replace_uniform'; + dm.params.mutation_scale=0.15; + dm.params.mutation_rate=0.08; + dm.params.replacement_type=''; + dm.params.below_limit=6; + dm.params.shrinkage_percentage=0.9; + dm.params.crossover_type='shuffle_random'; + dm.params.multi_point_binary=[]; + dm.params.multi_point_parameterized_binary=[]; + dm.params.multi_point_real=[]; + dm.params.shuffle_random=[]; + dm.params.num_parents=2; + dm.params.num_offspring=2; + dm.params.crossover_rate=0.8; + dm.params.fitness_type='merit_function'; + dm.params.constraint_penalty=1.0; + dm.params.replacement_type=''; + dm.params.convergence_type=false; + dm.params.num_generations=10; + dm.params.percent_change=0.1; + + case {'nl2sol'} + dm.type ='lstsq'; + dm.variables={'cdv','csv'}; + dm.lcspec ={}; + dm.responses={'lst'}; + dm.ghspec ={'grad'}; + dm.params.max_iterations=false; + dm.params.max_function_evaluations=false; + dm.params.convergence_tolerance=false; + dm.params.output=false; + dm.params.scaling=false; + dm.params.function_precision=1.e-10; + dm.params.absolute_conv_tol=-1.; + dm.params.x_conv_tol=-1.; + dm.params.singular_conv_tol=-1.; + dm.params.singular_radius=-1.; + dm.params.false_conv_tol=-1.; + dm.params.initial_trust_radius=-1.; + dm.params.covariance=0; + dm.params.regression_diagnostics=false; + case {'nlssol_sqp'} + dm.type ='lstsq'; + dm.variables={'cdv','csv'}; + dm.lcspec ={'lic','lec'}; + dm.responses={'lst','nic','nec'}; + dm.ghspec ={'grad'}; + dm.params.max_iterations=false; + dm.params.max_function_evaluations=false; + dm.params.convergence_tolerance=false; + dm.params.constraint_tolerance=false; + dm.params.output=false; + dm.params.speculative=false; + dm.params.scaling=false; + dm.params.verify_level=-1; + dm.params.function_precision=1.e-10; + dm.params.linesearch_tolerance=0.9; + case {'optpp_g_newton'} + dm.type ='lstsq'; + dm.variables={'cdv','csv'}; + dm.lcspec ={'lic','lec'}; + dm.responses={'lst','nic','nec'}; + dm.ghspec ={'grad'}; + dm.params.max_iterations=false; + dm.params.max_function_evaluations=false; + dm.params.convergence_tolerance=false; + dm.params.output=false; + dm.params.speculative=false; + dm.params.scaling=false; + dm.params.value_based_line_search=false; + dm.params.gradient_based_line_search=false; + dm.params.trust_region=false; + dm.params.tr_pds=false; + dm.params.max_step=1000.; + dm.params.gradient_tolerance=1.e-4; + dm.params.merit_function='argaez_tapia'; + dm.params.central_path=dm.params.merit_function; + dm.params.steplength_to_boundary=0.99995; + dm.params.centering_parameter=0.2; + + case {'nond_sampling'} + dm.type ='nond'; + dm.variables={'nuv','csv'}; + dm.lcspec ={}; + dm.responses={'rf'}; + dm.ghspec ={}; + dm.params.seed=false; + dm.params.fixed_seed=false; + dm.params.samples=false; + dm.params.sample_type='lhs'; + dm.params.all_variables=false; + dm.params.variance_based_decomp=false; + dm.params.previous_samples=0; + case {'nond_local_reliability'} + dm.type ='nond'; + dm.variables={'nuv','csv'}; + dm.lcspec ={}; + dm.responses={'rf'}; + dm.ghspec ={'grad'}; + dm.params.max_iterations=false; + dm.params.convergence_tolerance=false; + dm.params.mpp_search=false; + dm.params.sqp=false; + dm.params.nip=false; + dm.params.integration='first_order'; + dm.params.refinement=false; + dm.params.samples=0; + dm.params.seed=false; + case {'nond_global_reliability'} + dm.type ='nond'; + dm.variables={'nuv','csv'}; + dm.lcspec ={}; + dm.responses={'rf'}; + dm.ghspec ={'grad'}; + dm.params.x_gaussian_process=false; + dm.params.u_gaussian_process=false; + dm.params.all_variables=false; + dm.params.seed=false; + case {'nond_polynomial_chaos'} + dm.type ='nond'; + dm.variables={'nuv','csv'}; + dm.lcspec ={}; + dm.responses={'rf'}; + dm.ghspec ={'grad'}; + dm.params.expansion_order=[]; + dm.params.expansion_terms=[]; + dm.params.quadrature_order=[]; + dm.params.sparse_grid_level=[]; + dm.params.expansion_samples=[]; + dm.params.incremental_lhs=false; + dm.params.collocation_points=[]; + dm.params.collocation_ratio=[]; + dm.params.reuse_samples=false; + dm.params.expansion_import_file=''; + dm.params.seed=false; + dm.params.fixed_seed=false; + dm.params.samples=0; + dm.params.sample_type='lhs'; + dm.params.all_variables=false; + case {'nond_stoch_collocation'} + dm.type ='nond'; + dm.variables={'nuv','csv'}; + dm.lcspec ={}; + dm.responses={'rf'}; + dm.ghspec ={'grad'}; + dm.params.quadrature_order=[]; + dm.params.sparse_grid_level=[]; + dm.params.seed=false; + dm.params.fixed_seed=false; + dm.params.samples=0; + dm.params.sample_type='lhs'; + dm.params.all_variables=false; + case {'nond_evidence'} + dm.type ='nond'; + dm.variables={'nuv','csv'}; + dm.lcspec ={}; + dm.responses={'rf'}; + dm.ghspec ={'grad'}; + dm.params.seed=false; + dm.params.samples=10000; + + case {'dace'} + dm.type ='dace'; + dm.variables={'cdv','csv'}; + dm.lcspec ={}; + dm.responses={'rf'}; + dm.ghspec ={}; + dm.params.grid=false; + dm.params.random=false; + dm.params.oas=false; + dm.params.lhs=false; + dm.params.oa_lhs=false; + dm.params.box_behnken=false; + dm.params.central_composite=false; + dm.params.seed=false; + dm.params.fixed_seed=false; + dm.params.samples=false; + dm.params.symbols=false; + dm.params.quality_metrics=false; + dm.params.variance_based_decomp=false; + case {'fsu_quasi_mc'} + dm.type ='dace'; + dm.variables={'cdv','csv'}; + dm.lcspec ={}; + dm.responses={'rf'}; + dm.ghspec ={}; + dm.params.halton=false; + dm.params.hammersley=false; + dm.params.samples=0; + dm.params.sequence_start=[0]; + dm.params.sequence_leap=[1]; + dm.params.prime_base=false; + dm.params.fixed_sequence=false; + dm.params.latinize=false; + dm.params.variance_based_decomp=false; + dm.params.quality_metrics=false; + case {'fsu_cvt'} + dm.type ='dace'; + dm.variables={'cdv','csv'}; + dm.lcspec ={}; + dm.responses={'rf'}; + dm.ghspec ={}; + dm.params.seed=false; + dm.params.fixed_seed=false; + dm.params.samples=0; + dm.params.num_trials=10000; + dm.params.trial_type='random'; + dm.params.latinize=false; + dm.params.variance_based_decomp=false; + dm.params.quality_metrics=false; + + case {'vector_parameter_study'} + dm.type ='param'; + dm.variables={'cdv','csv'}; + dm.lcspec ={}; + dm.responses={'rf'}; + dm.ghspec ={}; + dm.params.output=false; + dm.params.final_point=[]; + dm.params.step_length=[]; + dm.params.num_steps=[]; + dm.params.step_vector=[]; + dm.params.num_steps=[]; + case {'list_parameter_study'} + dm.type ='param'; + dm.variables={'cdv','csv'}; + dm.lcspec ={}; + dm.responses={'rf'}; + dm.ghspec ={}; + dm.params.output=false; + dm.params.list_of_points=[]; + case {'centered_parameter_study'} + dm.type ='param'; + dm.variables={'cdv','csv'}; + dm.lcspec ={}; + dm.responses={'rf'}; + dm.ghspec ={}; + dm.params.output=false; + dm.params.percent_delta=[]; + dm.params.deltas_per_variable=[]; + case {'multidim_parameter_study'} + dm.type ='param'; + dm.variables={'cdv','csv'}; + dm.lcspec ={}; + dm.responses={'rf'}; + dm.ghspec ={}; + dm.params.output=false; + dm.params.partitions=[]; + + otherwise + error('Unimplemented method: ''%s''.',dm.method); + end + + end + +% if more than one argument, issue warning + + otherwise + warning('dakota_method:extra_arg',... + 'Extra arguments for object of class ''%s''.',... + class(dm)); + end + + end + end +end Index: /issm/trunk/src/m/classes/@dakota_method/display.m =================================================================== --- /issm/trunk/src/m/classes/@dakota_method/display.m (revision 27) +++ /issm/trunk/src/m/classes/@dakota_method/display.m (revision 27) @@ -0,0 +1,39 @@ +% +% display for the dakota_method class. +% +% []=display(dm) +% +function []=display(dm) + +if ~isa(dm,'dakota_method') + error('Object ''%s'' is a ''%s'' class object, not ''%s''.',... + inputname(1),class(dm),'dakota_method'); +end + +% display the object + +for i=1:numel(dm) + disp(sprintf('\nclass ''%s'' object ''%s%s'' = \n',... + class(dm),inputname(1),string_dim(dm,i))); + disp(sprintf(' method: ''%s''' ,dm(i).method)); + disp(sprintf(' type: ''%s''' ,dm(i).type)); + disp(sprintf(' variables: %s' ,string_cell(dm(i).variables))); + disp(sprintf(' lcspec: %s' ,string_cell(dm(i).lcspec))); + disp(sprintf(' responses: %s' ,string_cell(dm(i).responses))); + disp(sprintf(' ghspec: %s\n' ,string_cell(dm(i).ghspec))); + +% display the parameters within the object + + fnames=fieldnames(dm(i).params); + maxlen=0; + for j=1:numel(fnames) + maxlen=max(maxlen,length(fnames{j})); + end + + for j=1:numel(fnames) + disp(sprintf([' params.%-' num2str(maxlen+1) 's: %s'],... + fnames{j},any2str(dm(i).params.(fnames{j})))); + end +end + +end Index: /issm/trunk/src/m/classes/@dakota_method/dmeth_params_merge.m =================================================================== --- /issm/trunk/src/m/classes/@dakota_method/dmeth_params_merge.m (revision 27) +++ /issm/trunk/src/m/classes/@dakota_method/dmeth_params_merge.m (revision 27) @@ -0,0 +1,27 @@ +% +% merge a structure of parameters into a dakota_method object. +% +% [dm]=dmeth_params_merge(dm,params) +% +function [dm]=dmeth_params_merge(dm,params) + +if ~isa(dm,'dakota_method') + error('Object ''%s'' is a ''%s'' class object, not ''%s''.',... + inputname(1),class(dm),'dakota_method'); +end + +% loop through each parameter field in the structure + +fnames=fieldnames(params); + +for i=1:numel(fnames) + if isfield(dm.params,fnames{i}) + dm.params.(fnames{i})=params.(fnames{i}); + else + warning('dmeth_params_merge:unknown_param',... + 'No parameter ''%s'' for dakota_method ''%s''.',... + fnames{i},dm.method); + end +end + +end Index: /issm/trunk/src/m/classes/@dakota_method/dmeth_params_set.m =================================================================== --- /issm/trunk/src/m/classes/@dakota_method/dmeth_params_set.m (revision 27) +++ /issm/trunk/src/m/classes/@dakota_method/dmeth_params_set.m (revision 27) @@ -0,0 +1,25 @@ +% +% set parameters of a dakota_method object. +% +% [dm]=dmeth_params_set(dm,varargin) +% +function [dm]=dmeth_params_set(dm,varargin) + +if ~isa(dm,'dakota_method') + error('Object ''%s'' is a ''%s'' class object, not ''%s''.',... + inputname(1),class(dm),'dakota_method'); +end + +% loop through each parameter field in the input list + +for i=1:2:length(varargin) + if isfield(dm.params,varargin{i}) + dm.params.(varargin{i})=varargin{i+1}; + else + warning('dmeth_params_set:unknown_param',... + 'No parameter ''%s'' for dakota_method ''%s''.',... + varargin{i},dm.method); + end +end + +end Index: /issm/trunk/src/m/classes/@dakota_method/dmeth_params_write.m =================================================================== --- /issm/trunk/src/m/classes/@dakota_method/dmeth_params_write.m (revision 27) +++ /issm/trunk/src/m/classes/@dakota_method/dmeth_params_write.m (revision 27) @@ -0,0 +1,565 @@ +% +% write the parameters from a dakota_method object. +% +% []=dmeth_params_write(dm,fid,sbeg) +% +function []=dmeth_params_write(dm,fid,sbeg) + +if ~isa(dm,'dakota_method') + error('Object ''%s'' is a ''%s'' class object, not ''%s''.',... + inputname(1),class(dm),'dakota_method'); +end + +if ~exist('sbeg','var') + sbeg='\t '; +end + +% perform some error checking, but leave the rest to dakota. +% unfortunately this prevents merely looping through the fields +% of the parameters structure. + +% write method-independent controls + +% param_write(fid,sbeg,'id_method',' = ','\n',dm.params); +% param_write(fid,sbeg,'model_pointer',' = ','\n',dm.params); + +% write method-dependent controls + +switch dm.type + case {'dot'} + param_write(fid,sbeg,'max_iterations',' = ','\n',dm.params); + param_write(fid,sbeg,'max_function_evaluations',' = ','\n',dm.params); + param_write(fid,sbeg,'convergence_tolerance',' = ','\n',dm.params); + param_write(fid,sbeg,'constraint_tolerance',' = ','\n',dm.params); + param_write(fid,sbeg,'output',' ','\n',dm.params); + param_write(fid,sbeg,'speculative','','\n',dm.params); + param_write(fid,sbeg,'scaling','','\n',dm.params); + switch dm.method + case{'dot_bfgs',... + 'dot_frcg',... + 'dot_mmfd',... + 'dot_slp',... + 'dot_sqp'} + param_write(fid,sbeg,'optimization_type',' = ','\n',dm.params); + + otherwise + error('Unrecognized ''%s'' method: ''%s''.',dm.type,dm.method); + end + + case {'npsol'} + param_write(fid,sbeg,'max_iterations',' = ','\n',dm.params); + param_write(fid,sbeg,'max_function_evaluations',' = ','\n',dm.params); + param_write(fid,sbeg,'convergence_tolerance',' = ','\n',dm.params); + param_write(fid,sbeg,'constraint_tolerance',' = ','\n',dm.params); + param_write(fid,sbeg,'output',' ','\n',dm.params); + param_write(fid,sbeg,'speculative','','\n',dm.params); + param_write(fid,sbeg,'scaling','','\n',dm.params); + switch dm.method + case {'npsol_sqp'} + param_write(fid,sbeg,'verify_level',' = ','\n',dm.params); + param_write(fid,sbeg,'function_precision',' = ','\n',dm.params); + param_write(fid,sbeg,'linesearch_tolerance',' = ','\n',dm.params); + + otherwise + error('Unrecognized ''%s'' method: ''%s''.',dm.type,dm.method); + end + + case {'conmin'} + param_write(fid,sbeg,'max_iterations',' = ','\n',dm.params); + param_write(fid,sbeg,'max_function_evaluations',' = ','\n',dm.params); + param_write(fid,sbeg,'convergence_tolerance',' = ','\n',dm.params); + param_write(fid,sbeg,'constraint_tolerance',' = ','\n',dm.params); + param_write(fid,sbeg,'output',' ','\n',dm.params); + param_write(fid,sbeg,'speculative','','\n',dm.params); + param_write(fid,sbeg,'scaling','','\n',dm.params); + switch dm.method + case {'conmin_frcg',... + 'conmin_mfd'} + + otherwise + error('Unrecognized ''%s'' method: ''%s''.',dm.type,dm.method); + end + + case {'optpp'} + param_write(fid,sbeg,'max_iterations',' = ','\n',dm.params); + param_write(fid,sbeg,'max_function_evaluations',' = ','\n',dm.params); + param_write(fid,sbeg,'convergence_tolerance',' = ','\n',dm.params); + param_write(fid,sbeg,'output',' ','\n',dm.params); + param_write(fid,sbeg,'speculative','','\n',dm.params); + param_write(fid,sbeg,'scaling','','\n',dm.params); + switch dm.method + case {'optpp_cg'} + param_write(fid,sbeg,'max_step',' = ','\n',dm.params); + param_write(fid,sbeg,'gradient_tolerance',' = ','\n',dm.params); + + case {'optpp_q_newton',... + 'optpp_fd_newton',... + 'optpp_newton'} + if (dm.params.value_based_line_search + ... + dm.params.gradient_based_line_search + ... + dm.params.trust_region + ... + dm.params.tr_pds > 1) + error('''%s'' method must have only one algorithm.',... + dm.method); + end + param_write(fid,sbeg,'value_based_line_search','','\n',dm.params); + param_write(fid,sbeg,'gradient_based_line_search','','\n',dm.params); + param_write(fid,sbeg,'trust_region','','\n',dm.params); + param_write(fid,sbeg,'tr_pds','','\n',dm.params); + param_write(fid,sbeg,'max_step',' = ','\n',dm.params); + param_write(fid,sbeg,'gradient_tolerance',' = ','\n',dm.params); + param_write(fid,sbeg,'merit_function',' = ','\n',dm.params); + param_write(fid,sbeg,'central_path',' = ','\n',dm.params); + param_write(fid,sbeg,'steplength_to_boundary',' = ','\n',dm.params); + param_write(fid,sbeg,'centering_parameter',' = ','\n',dm.params); + + case {'optpp_pds'} + param_write(fid,sbeg,'search_scheme_size',' = ','\n',dm.params); + + otherwise + error('Unrecognized ''%s'' method: ''%s''.',dm.type,dm.method); + end + + case {'apps'} + param_write(fid,sbeg,'max_function_evaluations',' = ','\n',dm.params); + param_write(fid,sbeg,'constraint_tolerance',' = ','\n',dm.params); + param_write(fid,sbeg,'output',' ','\n',dm.params); + param_write(fid,sbeg,'scaling','','\n',dm.params); + switch dm.method + case {'asynch_pattern_search'} + param_write(fid,sbeg,'initial_delta',' = ','\n',dm.params); + param_write(fid,sbeg,'threshold_delta',' = ','\n',dm.params); + param_write(fid,sbeg,'contraction_factor',' = ','\n',dm.params); + param_write(fid,sbeg,'solution_target',' = ','\n',dm.params); + param_write(fid,sbeg,'synchronization',' = ','\n',dm.params); + param_write(fid,sbeg,'merit_function',' = ','\n',dm.params); + param_write(fid,sbeg,'constraint_penalty',' = ','\n',dm.params); + param_write(fid,sbeg,'smoothing_factor',' = ','\n',dm.params); + + otherwise + error('Unrecognized ''%s'' method: ''%s''.',dm.type,dm.method); + end + + case {'coliny'} + param_write(fid,sbeg,'max_iterations',' = ','\n',dm.params); + param_write(fid,sbeg,'max_function_evaluations',' = ','\n',dm.params); + param_write(fid,sbeg,'convergence_tolerance',' = ','\n',dm.params); + param_write(fid,sbeg,'output',' ','\n',dm.params); + param_write(fid,sbeg,'scaling','','\n',dm.params); + + param_write(fid,sbeg,'show_misc_options','','\n',dm.params); + param_write(fid,sbeg,'misc_options',' = ','\n',dm.params); + param_write(fid,sbeg,'solution_accuracy',' = ','\n',dm.params); + switch dm.method + case {'coliny_cobyla'} + param_write(fid,sbeg,'initial_delta',' = ','\n',dm.params); + param_write(fid,sbeg,'threshold_delta',' = ','\n',dm.params); + + case {'coliny_direct'} + param_write(fid,sbeg,'division',' = ','\n',dm.params); + param_write(fid,sbeg,'global_balance_parameter',' = ','\n',dm.params); + param_write(fid,sbeg,'local_balance_parameter',' = ','\n',dm.params); + param_write(fid,sbeg,'max_boxsize_limit',' = ','\n',dm.params); + param_write(fid,sbeg,'min_boxsize_limit',' = ','\n',dm.params); + param_write(fid,sbeg,'constraint_penalty',' = ','\n',dm.params); + + case {'coliny_ea'} + param_write(fid,sbeg,'seed',' = ','\n',dm.params); + param_write(fid,sbeg,'population_size',' = ','\n',dm.params); + param_write(fid,sbeg,'initialization_type',' = ','\n',dm.params); + param_write(fid,sbeg,'fitness_type',' = ','\n',dm.params); + param_write(fid,sbeg,'replacement_type',' = ','\n',dm.params); + param_write(fid,sbeg,'random',' = ','\n',dm.params); + param_write(fid,sbeg,'chc',' = ','\n',dm.params); + param_write(fid,sbeg,'elitist',' = ','\n',dm.params); + param_write(fid,sbeg,'new_solutions_generated',' = ','\n',dm.params); + param_write(fid,sbeg,'crossover_type',' = ','\n',dm.params); + param_write(fid,sbeg,'crossover_rate',' = ','\n',dm.params); + param_write(fid,sbeg,'mutation_type',' = ','\n',dm.params); + param_write(fid,sbeg,'mutation_scale',' = ','\n',dm.params); + param_write(fid,sbeg,'mutation_range',' = ','\n',dm.params); + param_write(fid,sbeg,'dimension_ratio',' = ','\n',dm.params); + param_write(fid,sbeg,'mutation_rate',' = ','\n',dm.params); + param_write(fid,sbeg,'non_adaptive','','\n',dm.params); + + case {'coliny_pattern_search'} + param_write(fid,sbeg,'stochastic','','\n',dm.params); + param_write(fid,sbeg,'seed',' = ','\n',dm.params); + param_write(fid,sbeg,'initial_delta',' = ','\n',dm.params); + param_write(fid,sbeg,'threshold_delta',' = ','\n',dm.params); + param_write(fid,sbeg,'constraint_penalty',' = ','\n',dm.params); + param_write(fid,sbeg,'constant_penalty','','\n',dm.params); + param_write(fid,sbeg,'pattern_basis',' = ','\n',dm.params); + param_write(fid,sbeg,'total_pattern_size',' = ','\n',dm.params); + param_write(fid,sbeg,'no_expansion','','\n',dm.params); + param_write(fid,sbeg,'expand_after_success',' = ','\n',dm.params); + param_write(fid,sbeg,'contraction_factor',' = ','\n',dm.params); + param_write(fid,sbeg,'synchronization',' = ','\n',dm.params); + param_write(fid,sbeg,'exploratory_moves',' = ','\n',dm.params); + + case {'coliny_solis_wets'} + param_write(fid,sbeg,'seed',' = ','\n',dm.params); + param_write(fid,sbeg,'initial_delta',' = ','\n',dm.params); + param_write(fid,sbeg,'threshold_delta',' = ','\n',dm.params); + param_write(fid,sbeg,'no_expansion','','\n',dm.params); + param_write(fid,sbeg,'expand_after_success',' = ','\n',dm.params); + param_write(fid,sbeg,'contract_after_failure',' = ','\n',dm.params); + param_write(fid,sbeg,'contraction_factor',' = ','\n',dm.params); + param_write(fid,sbeg,'constraint_penalty',' = ','\n',dm.params); + param_write(fid,sbeg,'constant_penalty','','\n',dm.params); + + otherwise + error('Unrecognized ''%s'' method: ''%s''.',dm.type,dm.method); + end + + case {'ncsu'} + param_write(fid,sbeg,'max_iterations',' = ','\n',dm.params); + param_write(fid,sbeg,'max_function_evaluations',' = ','\n',dm.params); + param_write(fid,sbeg,'scaling','','\n',dm.params); + switch dm.method + case {'ncsu_direct'} + param_write(fid,sbeg,'solution_accuracy',' = ','\n',dm.params); + param_write(fid,sbeg,'min_boxsize_limit',' = ','\n',dm.params); + param_write(fid,sbeg,'vol_boxsize_limit',' = ','\n',dm.params); + + otherwise + error('Unrecognized ''%s'' method: ''%s''.',dm.type,dm.method); + end + + case {'jega'} + param_write(fid,sbeg,'max_iterations',' = ','\n',dm.params); + param_write(fid,sbeg,'max_function_evaluations',' = ','\n',dm.params); + param_write(fid,sbeg,'output',' ','\n',dm.params); + param_write(fid,sbeg,'scaling','','\n',dm.params); + + param_write(fid,sbeg,'seed',' = ','\n',dm.params); + param_write(fid,sbeg,'log_file',' = ','\n',dm.params); + param_write(fid,sbeg,'population_size',' = ','\n',dm.params); + param_write(fid,sbeg,'print_each_pop','','\n',dm.params); + param_write(fid,sbeg,'output',' = ','\n',dm.params); + param_write(fid,sbeg,'initialization_type',' = ','\n',dm.params); + param_write(fid,sbeg,'mutation_type',' = ','\n',dm.params); + param_write(fid,sbeg,'mutation_scale',' = ','\n',dm.params); + param_write(fid,sbeg,'mutation_rate',' = ','\n',dm.params); + param_write(fid,sbeg,'replacement_type',' = ','\n',dm.params); + param_write(fid,sbeg,'below_limit',' = ','\n',dm.params); + param_write(fid,sbeg,'shrinkage_percentage',' = ','\n',dm.params); + param_write(fid,sbeg,'crossover_type',' = ','\n',dm.params); + param_write(fid,sbeg,'multi_point_binary',' = ','\n',dm.params); + param_write(fid,sbeg,'multi_point_parameterized_binary',' = ','\n',dm.params); + param_write(fid,sbeg,'multi_point_real',' = ','\n',dm.params); + param_write(fid,sbeg,'shuffle_random',' = ','\n',dm.params); + param_write(fid,sbeg,'num_parents',' = ','\n',dm.params); + param_write(fid,sbeg,'num_offspring',' = ','\n',dm.params); + param_write(fid,sbeg,'crossover_rate',' = ','\n',dm.params); + + switch dm.method + case {'moga'} + param_write(fid,sbeg,'fitness_type',' = ','\n',dm.params); + param_write(fid,sbeg,'niching_type',' = ','\n',dm.params); + if ~isempty(dm.params.radial) && ... + ~isempty(dm.params.distance) + error('''%s'' method must have only one niching distance.',... + dm.method); + end + param_write(fid,sbeg,'radial',' = ','\n',dm.params); + param_write(fid,sbeg,'distance',' = ','\n',dm.params); + param_write(fid,sbeg,'metric_tracker','','\n',dm.params); + param_write(fid,sbeg,'percent_change',' = ','\n',dm.params); + param_write(fid,sbeg,'num_generations',' = ','\n',dm.params); + param_write(fid,sbeg,'postprocessor_type',' = ','\n',dm.params); + param_write(fid,sbeg,'orthogonal_distance',' = ','\n',dm.params); + + case {'soga'} + param_write(fid,sbeg,'fitness_type',' = ','\n',dm.params); + param_write(fid,sbeg,'constraint_penalty',' = ','\n',dm.params); + param_write(fid,sbeg,'replacement_type',' = ','\n',dm.params); + param_write(fid,sbeg,'convergence_type',' = ','\n',dm.params); + param_write(fid,sbeg,'num_generations',' = ','\n',dm.params); + param_write(fid,sbeg,'percent_change',' = ','\n',dm.params); + + otherwise + error('Unrecognized ''%s'' method: ''%s''.',dm.type,dm.method); + end + + case {'lstsq'} + switch dm.method + case {'nl2sol'} + param_write(fid,sbeg,'max_iterations',' = ','\n',dm.params); + param_write(fid,sbeg,'max_function_evaluations',' = ','\n',dm.params); + param_write(fid,sbeg,'convergence_tolerance',' = ','\n',dm.params); + param_write(fid,sbeg,'output',' ','\n',dm.params); + param_write(fid,sbeg,'scaling','','\n',dm.params); + + param_write(fid,sbeg,'function_precision',' = ','\n',dm.params); + param_write(fid,sbeg,'absolute_conv_tol',' = ','\n',dm.params); + param_write(fid,sbeg,'x_conv_tol',' = ','\n',dm.params); + param_write(fid,sbeg,'singular_conv_tol',' = ','\n',dm.params); + param_write(fid,sbeg,'singular_radius',' = ','\n',dm.params); + param_write(fid,sbeg,'false_conv_tol',' = ','\n',dm.params); + param_write(fid,sbeg,'initial_trust_radius',' = ','\n',dm.params); + param_write(fid,sbeg,'covariance',' = ','\n',dm.params); + param_write(fid,sbeg,'regression_diagnostics','','\n',dm.params); + + case {'nlssol_sqp'} + param_write(fid,sbeg,'max_iterations',' = ','\n',dm.params); + param_write(fid,sbeg,'max_function_evaluations',' = ','\n',dm.params); + param_write(fid,sbeg,'convergence_tolerance',' = ','\n',dm.params); + param_write(fid,sbeg,'constraint_tolerance',' = ','\n',dm.params); + param_write(fid,sbeg,'output',' ','\n',dm.params); + param_write(fid,sbeg,'speculative','','\n',dm.params); + param_write(fid,sbeg,'scaling','','\n',dm.params); + + param_write(fid,sbeg,'verify_level',' = ','\n',dm.params); + param_write(fid,sbeg,'function_precision',' = ','\n',dm.params); + param_write(fid,sbeg,'linesearch_tolerance',' = ','\n',dm.params); + + case {'optpp_g_newton'} + param_write(fid,sbeg,'max_iterations',' = ','\n',dm.params); + param_write(fid,sbeg,'max_function_evaluations',' = ','\n',dm.params); + param_write(fid,sbeg,'convergence_tolerance',' = ','\n',dm.params); + param_write(fid,sbeg,'output',' ','\n',dm.params); + param_write(fid,sbeg,'speculative','','\n',dm.params); + param_write(fid,sbeg,'scaling','','\n',dm.params); + + if (dm.params.value_based_line_search + ... + dm.params.gradient_based_line_search + ... + dm.params.trust_region + ... + dm.params.tr_pds > 1) + error('''%s'' method must have only one algorithm.',... + dm.method); + end + param_write(fid,sbeg,'value_based_line_search','','\n',dm.params); + param_write(fid,sbeg,'gradient_based_line_search','','\n',dm.params); + param_write(fid,sbeg,'trust_region','','\n',dm.params); + param_write(fid,sbeg,'tr_pds','','\n',dm.params); + param_write(fid,sbeg,'max_step',' = ','\n',dm.params); + param_write(fid,sbeg,'gradient_tolerance',' = ','\n',dm.params); + param_write(fid,sbeg,'merit_function',' = ','\n',dm.params); + param_write(fid,sbeg,'central_path',' = ','\n',dm.params); + param_write(fid,sbeg,'steplength_to_boundary',' = ','\n',dm.params); + param_write(fid,sbeg,'centering_parameter',' = ','\n',dm.params); + + otherwise + error('Unrecognized ''%s'' method: ''%s''.',dm.type,dm.method); + end + + case {'nond'} + switch dm.method + case {'nond_sampling'} + param_write(fid,sbeg,'seed',' = ','\n',dm.params); + param_write(fid,sbeg,'fixed_seed','','\n',dm.params); + param_write(fid,sbeg,'samples',' = ','\n',dm.params); + param_write(fid,sbeg,'sample_type',' = ','\n',dm.params); + param_write(fid,sbeg,'all_variables','','\n',dm.params); + param_write(fid,sbeg,'variance_based_decomp','','\n',dm.params); + param_write(fid,sbeg,'previous_samples',' = ','\n',dm.params); + + case {'nond_local_reliability'} + param_write(fid,sbeg,'max_iterations',' = ','\n',dm.params); + param_write(fid,sbeg,'convergence_tolerance',' = ','\n',dm.params); + + param_write(fid,sbeg,'mpp_search',' = ','\n',dm.params); + if (dm.params.sqp + ... + dm.params.nip > 1) + error('''%s'' method must have only one algorithm.',... + dm.method); + end + param_write(fid,sbeg,'sqp','','\n',dm.params); + param_write(fid,sbeg,'nip','','\n',dm.params); + param_write(fid,sbeg,'integration',' = ','\n',dm.params); + param_write(fid,sbeg,'refinement',' = ','\n',dm.params); + param_write(fid,sbeg,'samples',' = ','\n',dm.params); + param_write(fid,sbeg,'seed',' = ','\n',dm.params); + + case {'nond_global_reliability'} + if (dm.params.x_gaussian_process + ... + dm.params.u_gaussian_process ~= 1) + error('''%s'' method must have one and only one algorithm.',... + dm.method); + end + param_write(fid,sbeg,'x_gaussian_process','','\n',dm.params); + param_write(fid,sbeg,'u_gaussian_process','','\n',dm.params); + param_write(fid,sbeg,'all_variables','','\n',dm.params); + param_write(fid,sbeg,'seed',' = ','\n',dm.params); + + case {'nond_polynomial_chaos'} + param_write(fid,sbeg,'expansion_order',' = ','\n',dm.params); + param_write(fid,sbeg,'expansion_terms',' = ','\n',dm.params); + param_write(fid,sbeg,'quadrature_order',' = ','\n',dm.params); + param_write(fid,sbeg,'sparse_grid_level',' = ','\n',dm.params); + param_write(fid,sbeg,'expansion_samples',' = ','\n',dm.params); + param_write(fid,sbeg,'incremental_lhs','','\n',dm.params); + param_write(fid,sbeg,'collocation_points',' = ','\n',dm.params); + param_write(fid,sbeg,'collocation_ratio',' = ','\n',dm.params); + param_write(fid,sbeg,'reuse_samples','','\n',dm.params); + param_write(fid,sbeg,'expansion_import_file',' = ','\n',dm.params); + param_write(fid,sbeg,'seed',' = ','\n',dm.params); + param_write(fid,sbeg,'fixed_seed','','\n',dm.params); + param_write(fid,sbeg,'samples',' = ','\n',dm.params); + param_write(fid,sbeg,'sample_type',' = ','\n',dm.params); + param_write(fid,sbeg,'all_variables','','\n',dm.params); + + case {'nond_stoch_collocation'} + param_write(fid,sbeg,'quadrature_order',' = ','\n',dm.params); + param_write(fid,sbeg,'sparse_grid_level',' = ','\n',dm.params); + param_write(fid,sbeg,'seed',' = ','\n',dm.params); + param_write(fid,sbeg,'fixed_seed','','\n',dm.params); + param_write(fid,sbeg,'samples',' = ','\n',dm.params); + param_write(fid,sbeg,'sample_type',' = ','\n',dm.params); + param_write(fid,sbeg,'all_variables','','\n',dm.params); + + case {'nond_evidence'} + param_write(fid,sbeg,'seed',' = ','\n',dm.params); + param_write(fid,sbeg,'samples',' = ','\n',dm.params); + + otherwise + error('Unrecognized ''%s'' method: ''%s''.',dm.type,dm.method); + end + + case {'dace'} + switch dm.method + case {'dace'} + if (dm.params.grid + ... + dm.params.random + ... + dm.params.oas + ... + dm.params.lhs + ... + dm.params.oa_lhs + ... + dm.params.box_behnken + ... + dm.params.central_composite ~= 1) + error('''%s'' method must have one and only one algorithm.',... + dm.method); + end + param_write(fid,sbeg,'grid','','\n',dm.params); + param_write(fid,sbeg,'random','','\n',dm.params); + param_write(fid,sbeg,'oas','','\n',dm.params); + param_write(fid,sbeg,'lhs','','\n',dm.params); + param_write(fid,sbeg,'oa_lhs','','\n',dm.params); + param_write(fid,sbeg,'box_behnken','','\n',dm.params); + param_write(fid,sbeg,'central_composite','','\n',dm.params); + param_write(fid,sbeg,'seed',' = ','\n',dm.params); + param_write(fid,sbeg,'fixed_seed','','\n',dm.params); + param_write(fid,sbeg,'samples',' = ','\n',dm.params); + param_write(fid,sbeg,'symbols',' = ','\n',dm.params); + param_write(fid,sbeg,'quality_metrics','','\n',dm.params); + param_write(fid,sbeg,'variance_based_decomp','','\n',dm.params); + + case {'fsu_quasi_mc'} + if (dm.params.halton + ... + dm.params.hammersley ~= 1) + error('''%s'' method must have one and only one sequence type.',... + dm.method); + end + param_write(fid,sbeg,'halton','','\n',dm.params); + param_write(fid,sbeg,'hammersley','','\n',dm.params); + param_write(fid,sbeg,'samples',' = ','\n',dm.params); + param_write(fid,sbeg,'sequence_start',' = ','\n',dm.params); + param_write(fid,sbeg,'sequence_leap',' = ','\n',dm.params); + param_write(fid,sbeg,'prime_base',' = ','\n',dm.params); + param_write(fid,sbeg,'fixed_sequence','','\n',dm.params); + param_write(fid,sbeg,'latinize','','\n',dm.params); + param_write(fid,sbeg,'variance_based_decomp','','\n',dm.params); + param_write(fid,sbeg,'quality_metrics','','\n',dm.params); + + case {'fsu_cvt'} + param_write(fid,sbeg,'seed',' = ','\n',dm.params); + param_write(fid,sbeg,'fixed_seed','','\n',dm.params); + param_write(fid,sbeg,'samples',' = ','\n',dm.params); + param_write(fid,sbeg,'num_trials',' = ','\n',dm.params); + param_write(fid,sbeg,'trial_type',' = ','\n',dm.params); + param_write(fid,sbeg,'latinize','','\n',dm.params); + param_write(fid,sbeg,'variance_based_decomp','','\n',dm.params); + param_write(fid,sbeg,'quality_metrics','','\n',dm.params); + + otherwise + error('Unrecognized ''%s'' method: ''%s''.',dm.type,dm.method); + end + + case {'param'} + param_write(fid,sbeg,'output',' ','\n',dm.params); + switch dm.method + case {'vector_parameter_study'} + if ~xor(isempty(dm.params.final_point), ... + isempty(dm.params.step_vector)) + error('''%s'' method must have one and only one specification.',... + dm.method); + end + if ~isempty(dm.params.final_point) + param_write(fid,sbeg,'final_point',' = ','\n',dm.params); + param_write(fid,sbeg,'step_length',' = ','\n',dm.params); + param_write(fid,sbeg,'num_steps',' = ','\n',dm.params); + elseif ~isempty(dm.params.step_vector) + param_write(fid,sbeg,'step_vector',' = ','\n',dm.params); + param_write(fid,sbeg,'num_steps',' = ','\n',dm.params); + end + + case {'list_parameter_study'} + param_write(fid,sbeg,'list_of_points',' = ','\n',dm.params); + + case {'centered_parameter_study'} + param_write(fid,sbeg,'percent_delta',' = ','\n',dm.params); + param_write(fid,sbeg,'deltas_per_variable',' = ','\n',dm.params); + + case {'multidim_parameter_study'} + param_write(fid,sbeg,'partitions',' = ','\n',dm.params); + + otherwise + error('Unrecognized ''%s'' method: ''%s''.',dm.type,dm.method); + end + + otherwise + error('Unrecognized method type: ''%s''.',dm.type); +end + +end + +%% function to write a structure of parameters + +function []=param_struc_write(fidi,sbeg,smid,send,params) + +% loop through each parameter field in the structure + +fnames=fieldnames(params); + +for i=1:numel(fnames) + param_write(fidi,sbeg,fnames{i},smid,send,params); +end + +end + +%% function to write a parameter + +function []=param_write(fidi,sbeg,pname,smid,send,params) + +% check for errors + +if ~isfield(params,pname) + warning('param_write:param_not_found',... + 'Parameter ''%s'' not found in ''%s''.',... + pname,inputname(6)); + return +elseif islogical(params.(pname)) && ~params.(pname) + return +elseif isempty(params.(pname)) + warning('param_write:param_empty',... + 'Parameter ''%s'' requires input of type ''%s''.',... + pname,class(params.(pname))); + return +end + +% construct the parameter string based on type + +if islogical(params.(pname)) + fprintf(fidi,[sbeg '%s' send],pname); +elseif ischar (params.(pname)) + fprintf(fidi,[sbeg '%s' smid '%s' send],pname,params.(pname)); +elseif isnumeric(params.(pname)) + fprintf(fidi,[sbeg '%s' smid '%g' send],pname,params.(pname)); +else + warning('param_write:param_unrecog',... + 'Parameter ''%s'' is of unrecognized type ''%s''.',... + pname,class(params.(pname))); + return +end + +end Index: /issm/trunk/src/m/classes/@least_squares_term/least_squares_term.m =================================================================== --- /issm/trunk/src/m/classes/@least_squares_term/least_squares_term.m (revision 26) +++ /issm/trunk/src/m/classes/@least_squares_term/least_squares_term.m (revision 27) @@ -48,35 +48,39 @@ end - function [desc] =dresp_desc(lst) + function [desc] =prop_desc(lst) desc=cell(size(lst)); for i=1:numel(lst) desc(i)=cellstr(lst(i).descriptor); end + desc=allempty(desc); end - function [stype ]=dresp_stype(lst) + function [stype] =prop_stype(lst) stype=cell(size(lst)); for i=1:numel(lst) stype(i)=cellstr(lst(i).scale_type); end + stype=allequal(stype,'none'); end - function [scale] =dresp_scale(lst) + function [scale] =prop_scale(lst) scale=zeros(size(lst)); for i=1:numel(lst) scale(i)=lst(i).scale; end + scale=allequal(scale,1.); end - function [weight]=dresp_weight(lst) + function [weight]=prop_weight(lst) weight=zeros(size(lst)); for i=1:numel(lst) weight(i)=lst(i).weight; end + weight=allequal(weight,1.); end - function [lower] =dresp_lower(lst) + function [lower] =prop_lower(lst) lower=[]; end - function [upper] =dresp_upper(lst) + function [upper] =prop_upper(lst) upper=[]; end - function [target]=dresp_target(lst) + function [target]=prop_target(lst) target=[]; end Index: /issm/trunk/src/m/classes/@linear_equality_constraint/display.m =================================================================== --- /issm/trunk/src/m/classes/@linear_equality_constraint/display.m (revision 27) +++ /issm/trunk/src/m/classes/@linear_equality_constraint/display.m (revision 27) @@ -0,0 +1,26 @@ +% +% display for the linear_equality_constraint class. +% +% []=display(lec) +% +function []=display(lec) + +if ~isa(lec,'linear_equality_constraint') + error('Object ''%s'' is a ''%s'' class object, not ''%s''.',... + inputname(1),class(lec),'linear_equality_constraint'); +end + +% display the object + +disp(sprintf('\n')); +for i=1:numel(lec) + disp(sprintf('class ''%s'' object ''%s%s'' = \n',... + class(lec),inputname(1),string_dim(lec,i))); + disp(sprintf(' matrix: %s' ,string_vec(lec(i).matrix))); + disp(sprintf(' target: %g' ,lec(i).target)); + disp(sprintf(' scale_type: ''%s''' ,lec(i).scale_type)); + disp(sprintf(' scale: %g\n' ,lec(i).scale)); +end + +end + Index: /issm/trunk/src/m/classes/@linear_equality_constraint/linear_equality_constraint.m =================================================================== --- /issm/trunk/src/m/classes/@linear_equality_constraint/linear_equality_constraint.m (revision 27) +++ /issm/trunk/src/m/classes/@linear_equality_constraint/linear_equality_constraint.m (revision 27) @@ -0,0 +1,95 @@ +% +% constructor for the linear_equality_constraint class. +% +% [lec]=linear_equality_constraint(varargin) +% +classdef linear_equality_constraint + properties + matrix = NaN; + target = 0.; + scale_type='none'; + scale = 1.; + end + + methods + function [lec]=linear_equality_constraint(varargin) + + switch nargin + +% create a default object + + case 0 + +% copy the object + + case 1 + if isa(varargin{1},'linear_equality_constraint') + lec=varargin{1}; + else + error('Object ''%s'' is a ''%s'' class object, not ''%s''.',... + inputname(1),class(varargin{1}),'linear_equality_constraint'); + end + +% create the object from the input + + otherwise + if (size(varargin{1},1) > 1) + warning('linear_equality_constraint:matrix_rows',... + 'Matrix for object of class ''%s'' has %d rows.',... + class(lec),size(varargin{1},1)); + end + lec.matrix =varargin{1}(1,:); + + if (nargin >= 2) + lec.target =varargin{2}; + if (nargin >= 3) + lec.scale_type=varargin{3}; + if (nargin >= 4) + lec.scale =varargin{4}; + + if (nargin > 4) + warning('linear_equality_constraint:extra_arg',... + 'Extra arguments for object of class ''%s''.',... + class(lec)); + end + end + end + end + end + end + + function [matrix]=prop_matrix(lec) + matrix=zeros(numel(lec),0); + for i=1:numel(lec) + matrix(i,1:size(lec(i).matrix,2))=lec(i).matrix(1,:); + end + end + function [lower] =prop_lower(lec) + lower=[]; + end + function [upper] =prop_upper(lec) + upper=[]; + end + function [target]=prop_target(lec) + target=zeros(size(lec)); + for i=1:numel(lec) + target(i)=lec(i).target; + end + target=allequal(target,0.); + end + function [stype] =prop_stype(lec) + stype=cell(size(lec)); + for i=1:numel(lec) + stype(i)=cellstr(lec(i).scale_type); + end + stype=allequal(stype,'none'); + end + function [scale] =prop_scale(lec) + scale=zeros(size(lec)); + for i=1:numel(lec) + scale(i)=lec(i).scale; + end + scale=allequal(scale,1.); + end + end +end Index: /issm/trunk/src/m/classes/@linear_inequality_constraint/display.m =================================================================== --- /issm/trunk/src/m/classes/@linear_inequality_constraint/display.m (revision 27) +++ /issm/trunk/src/m/classes/@linear_inequality_constraint/display.m (revision 27) @@ -0,0 +1,27 @@ +% +% display for the linear_inequality_constraint class. +% +% []=display(lic) +% +function []=display(lic) + +if ~isa(lic,'linear_inequality_constraint') + error('Object ''%s'' is a ''%s'' class object, not ''%s''.',... + inputname(1),class(lic),'linear_inequality_constraint'); +end + +% display the object + +disp(sprintf('\n')); +for i=1:numel(lic) + disp(sprintf('class ''%s'' object ''%s%s'' = \n',... + class(lic),inputname(1),string_dim(lic,i))); + disp(sprintf(' matrix: %s' ,string_vec(lic(i).matrix))); + disp(sprintf(' lower: %g' ,lic(i).lower)); + disp(sprintf(' upper: %g' ,lic(i).upper)); + disp(sprintf(' scale_type: ''%s''' ,lic(i).scale_type)); + disp(sprintf(' scale: %g\n' ,lic(i).scale)); +end + +end + Index: /issm/trunk/src/m/classes/@linear_inequality_constraint/linear_inequality_constraint.m =================================================================== --- /issm/trunk/src/m/classes/@linear_inequality_constraint/linear_inequality_constraint.m (revision 27) +++ /issm/trunk/src/m/classes/@linear_inequality_constraint/linear_inequality_constraint.m (revision 27) @@ -0,0 +1,104 @@ +% +% constructor for the linear_inequality_constraint class. +% +% [lic]=linear_inequality_constraint(varargin) +% +classdef linear_inequality_constraint + properties + matrix = NaN; + lower =-Inf; + upper = 0.; + scale_type='none'; + scale = 1.; + end + + methods + function [lic]=linear_inequality_constraint(varargin) + + switch nargin + +% create a default object + + case 0 + +% copy the object + + case 1 + if isa(varargin{1},'linear_inequality_constraint') + lic=varargin{1}; + else + error('Object ''%s'' is a ''%s'' class object, not ''%s''.',... + inputname(1),class(varargin{1}),'linear_inequality_constraint'); + end + +% create the object from the input + + otherwise + if (size(varargin{1},1) > 1) + warning('linear_inequality_constraint:matrix_rows',... + 'Matrix for object of class ''%s'' has %d rows.',... + class(lic),size(varargin{1},1)); + end + lic.matrix =varargin{1}(1,:); + + if (nargin >= 2) + lic.lower =varargin{2}; + if (nargin >= 3) + lic.upper =varargin{3}; + if (nargin >= 4) + lic.scale_type=varargin{4}; + if (nargin >= 5) + lic.scale =varargin{5}; + + if (nargin > 5) + warning('linear_inequality_constraint:extra_arg',... + 'Extra arguments for object of class ''%s''.',... + class(lic)); + end + end + end + end + end + end + end + + function [matrix]=prop_matrix(lic) + matrix=zeros(numel(lic),0); + for i=1:numel(lic) + matrix(i,1:size(lic(i).matrix,2))=lic(i).matrix(1,:); + end + end + function [lower] =prop_lower(lic) + lower=zeros(size(lic)); + for i=1:numel(lic) + lower(i)=lic(i).lower; + end + lower=allequal(lower,-Inf); + end + function [upper] =prop_upper(lic) + upper=zeros(size(lic)); + for i=1:numel(lic) + upper(i)=lic(i).upper; + end + upper=allequal(upper,0.); + end + function [target]=prop_target(lic) + target=[]; + end + function [stype] =prop_stype(lic) + stype=cell(size(lic)); + for i=1:numel(lic) + stype(i)=cellstr(lic(i).scale_type); + end + stype=allequal(stype,'none'); + end + function [scale] =prop_scale(lic) + scale=zeros(size(lic)); + for i=1:numel(lic) + scale(i)=lic(i).scale; + end + scale=allequal(scale,1.); + end + end +end + Index: /issm/trunk/src/m/classes/@model/model.m =================================================================== --- /issm/trunk/src/m/classes/@model/model.m (revision 26) +++ /issm/trunk/src/m/classes/@model/model.m (revision 27) @@ -88,6 +88,6 @@ %Materials parameters - md.rho_ice=0; - md.rho_water=0; + md.rho_ice=917; + md.rho_water=1023; md.heatcapacity=2093; md.latentheat=3.34*10^5; %(J/kg); @@ -99,5 +99,5 @@ %Physical parameters - md.g=0; + md.g=9.81; md.yts=365*24*3600; md.drag=NaN; @@ -148,8 +148,8 @@ %Statics parameters - md.eps_rel=0; - md.eps_abs=0; + md.eps_rel=0.01; + md.eps_abs=10; md.acceleration=0; - md.sparsity=0; + md.sparsity=0.001; md.connectivity=10; md.lowmem=0; @@ -172,5 +172,5 @@ %Control - md.control_type='drag'; + md.control_type={'drag'}; md.cont_vx=NaN; md.cont_vy=NaN; @@ -181,5 +181,5 @@ md.nsteps=0; md.maxiter=[]; - md.tolx=0; + md.tolx=10^-4; md.optscal=[]; md.mincontrolconstraint=0; @@ -188,5 +188,5 @@ md.epsvel=eps; md.meanvel=0; - + %Output parameters md.parameteroutput={}; @@ -199,5 +199,5 @@ md.strainrate=NaN; md.plot=1; - + %debugging md.debug=1; @@ -247,6 +247,6 @@ %Cielo parameters md.solverstring=' -mat_type aijmumps -ksp_type preonly -pc_type lu -mat_mumps_icntl_14 40 '; - - %Needed for running in parallel + + %Trash, still to organize md.analysis_type=''; @@ -260,22 +260,14 @@ %qmu - md.variables=struct([]); - md.responses=struct([]); + md.variables =struct(); + md.responses =struct(); + md.qmu_method=struct(); + md.qmu_params=struct(); md.dakotaresults=NaN; - md.dakotain =NaN; - md.dakotaout=NaN; - md.dakotadat=NaN; -% md.method ='sampling'; - md.seed =1234; - md.samples =10; - md.sample_type ='lhs'; - md.method ='local reliability'; - md.interval_type ='forward'; - md.fd_gradient_step_size= 0.01; - md.evaluation_concurrency= 2; - md.npart = 1; - - md.analysis_driver =''; - md.analysis_components=''; + md.dakotain =NaN; + md.dakotaout =NaN; + md.dakotadat =NaN; + + md.npart =1; %Solver options Index: /issm/trunk/src/m/classes/@nonlinear_equality_constraint/nonlinear_equality_constraint.m =================================================================== --- /issm/trunk/src/m/classes/@nonlinear_equality_constraint/nonlinear_equality_constraint.m (revision 26) +++ /issm/trunk/src/m/classes/@nonlinear_equality_constraint/nonlinear_equality_constraint.m (revision 27) @@ -7,5 +7,5 @@ properties descriptor=''; - target = NaN; + target = 0.; scale_type='none'; scale = 1.; @@ -52,36 +52,40 @@ end - function [desc] =dresp_desc(nec) + function [desc] =prop_desc(nec) desc=cell(size(nec)); for i=1:numel(nec) desc(i)=cellstr(nec(i).descriptor); end + desc=allempty(desc); end - function [stype ]=dresp_stype(nec) + function [stype] =prop_stype(nec) stype=cell(size(nec)); for i=1:numel(nec) stype(i)=cellstr(nec(i).scale_type); end + stype=allequal(stype,'none'); end - function [scale] =dresp_scale(nec) + function [scale] =prop_scale(nec) scale=zeros(size(nec)); for i=1:numel(nec) scale(i)=nec(i).scale; end + scale=allequal(scale,1.); end - function [weight]=dresp_weight(nec) + function [weight]=prop_weight(nec) weight=[]; end - function [lower] =dresp_lower(nec) + function [lower] =prop_lower(nec) lower=[]; end - function [upper] =dresp_upper(nec) + function [upper] =prop_upper(nec) upper=[]; end - function [target]=dresp_target(nec) + function [target]=prop_target(nec) target=zeros(size(nec)); for i=1:numel(nec) target(i)=nec(i).target; end + target=allequal(target,0.); end end Index: /issm/trunk/src/m/classes/@nonlinear_inequality_constraint/nonlinear_inequality_constraint.m =================================================================== --- /issm/trunk/src/m/classes/@nonlinear_inequality_constraint/nonlinear_inequality_constraint.m (revision 26) +++ /issm/trunk/src/m/classes/@nonlinear_inequality_constraint/nonlinear_inequality_constraint.m (revision 27) @@ -7,6 +7,6 @@ properties descriptor=''; - lower = NaN; - upper = NaN; + lower =-Inf; + upper = 0.; scale_type='none'; scale = 1.; @@ -54,38 +54,43 @@ end - function [desc] =dresp_desc(nic) + function [desc] =prop_desc(nic) desc=cell(size(nic)); for i=1:numel(nic) desc(i)=cellstr(nic(i).descriptor); end + desc=allempty(desc); end - function [stype ]=dresp_stype(nic) + function [stype] =prop_stype(nic) stype=cell(size(nic)); for i=1:numel(nic) stype(i)=cellstr(nic(i).scale_type); end + stype=allequal(stype,'none'); end - function [scale] =dresp_scale(nic) + function [scale] =prop_scale(nic) scale=zeros(size(nic)); for i=1:numel(nic) scale(i)=nic(i).scale; end + scale=allequal(scale,1.); end - function [weight]=dresp_weight(nic) + function [weight]=prop_weight(nic) weight=[]; end - function [lower] =dresp_lower(nic) + function [lower] =prop_lower(nic) lower=zeros(size(nic)); for i=1:numel(nic) lower(i)=nic(i).lower; end + lower=allequal(lower,-Inf); end - function [upper] =dresp_upper(nic) + function [upper] =prop_upper(nic) upper=zeros(size(nic)); for i=1:numel(nic) upper(i)=nic(i).upper; end + upper=allequal(upper,0.); end - function [target]=dresp_target(nic) + function [target]=prop_target(nic) target=[]; end Index: /issm/trunk/src/m/classes/@normal_uncertain/normal_uncertain.m =================================================================== --- /issm/trunk/src/m/classes/@normal_uncertain/normal_uncertain.m (revision 26) +++ /issm/trunk/src/m/classes/@normal_uncertain/normal_uncertain.m (revision 27) @@ -9,6 +9,6 @@ mean = NaN; stddev = NaN; - lower = NaN; - upper = NaN; + lower =-Inf; + upper = Inf; end @@ -44,4 +44,5 @@ nuv.mean =varargin{2}; nuv.stddev =varargin{3}; + if (nargin >= 4) nuv.lower =varargin{4}; @@ -58,5 +59,5 @@ end - function [desc] =dvar_desc(nuv) + function [desc] =prop_desc(nuv) desc=cell(size(nuv)); for i=1:numel(nuv) @@ -65,22 +66,22 @@ desc=allempty(desc); end - function [initpt]=dvar_initpt(nuv) + function [initpt]=prop_initpt(nuv) initpt=[]; end - function [lower] =dvar_lower(nuv) + function [lower] =prop_lower(nuv) lower=zeros(size(nuv)); for i=1:numel(nuv) lower(i)=nuv(i).lower; end - lower=allnan(lower); + lower=allequal(lower,-Inf); end - function [upper] =dvar_upper(nuv) + function [upper] =prop_upper(nuv) upper=zeros(size(nuv)); for i=1:numel(nuv) upper(i)=nuv(i).upper; end - upper=allnan(upper); + upper=allequal(upper, Inf); end - function [mean] =dvar_mean(nuv) + function [mean] =prop_mean(nuv) mean=zeros(size(nuv)); for i=1:numel(nuv) @@ -88,5 +89,5 @@ end end - function [stddev]=dvar_stddev(nuv) + function [stddev]=prop_stddev(nuv) stddev=zeros(size(nuv)); for i=1:numel(nuv) @@ -94,6 +95,12 @@ end end - function [initst]=dvar_initst(nuv) + function [initst]=prop_initst(nuv) initst=[]; + end + function [stype] =prop_stype(nuv) + stype={}; + end + function [scale] =prop_scale(nuv) + scale=[]; end end Index: /issm/trunk/src/m/classes/@objective_function/objective_function.m =================================================================== --- /issm/trunk/src/m/classes/@objective_function/objective_function.m (revision 26) +++ /issm/trunk/src/m/classes/@objective_function/objective_function.m (revision 27) @@ -48,35 +48,39 @@ end - function [desc] =dresp_desc(of) + function [desc] =prop_desc(of) desc=cell(size(of)); for i=1:numel(of) desc(i)=cellstr(of(i).descriptor); end + desc=allempty(desc); end - function [stype ]=dresp_stype(of) + function [stype] =prop_stype(of) stype=cell(size(of)); for i=1:numel(of) stype(i)=cellstr(of(i).scale_type); end + stype=allequal(stype,'none'); end - function [scale] =dresp_scale(of) + function [scale] =prop_scale(of) scale=zeros(size(of)); for i=1:numel(of) scale(i)=of(i).scale; end + scale=allequal(scale,1.); end - function [weight]=dresp_weight(of) + function [weight]=prop_weight(of) weight=zeros(size(of)); for i=1:numel(of) weight(i)=of(i).weight; end + weight=allequal(weight,1.); end - function [lower] =dresp_lower(of) + function [lower] =prop_lower(of) lower=[]; end - function [upper] =dresp_upper(of) + function [upper] =prop_upper(of) upper=[]; end - function [target]=dresp_target(of) + function [target]=prop_target(of) target=[]; end Index: /issm/trunk/src/m/classes/@response_function/response_function.m =================================================================== --- /issm/trunk/src/m/classes/@response_function/response_function.m (revision 26) +++ /issm/trunk/src/m/classes/@response_function/response_function.m (revision 27) @@ -52,29 +52,30 @@ end - function [desc] =dresp_desc(rf) + function [desc] =prop_desc(rf) desc=cell(size(rf)); for i=1:numel(rf) desc(i)=cellstr(rf(i).descriptor); end + desc=allempty(desc); end - function [stype ]=dresp_stype(rf) + function [stype] =prop_stype(rf) stype={}; end - function [scale] =dresp_scale(rf) + function [scale] =prop_scale(rf) scale=[]; end - function [weight]=dresp_weight(rf) + function [weight]=prop_weight(rf) weight=[]; end - function [lower] =dresp_lower(rf) + function [lower] =prop_lower(rf) lower=[]; end - function [upper] =dresp_upper(rf) + function [upper] =prop_upper(rf) upper=[]; end - function [target]=dresp_target(rf) + function [target]=prop_target(rf) target=[]; end - function [respl,probl,rell,grell]=dresp_levels(rf) + function [respl,probl,rell,grell]=prop_levels(rf) respl=cell(size(rf)); probl=cell(size(rf)); Index: /issm/trunk/src/m/classes/public/SectionValues.m =================================================================== --- /issm/trunk/src/m/classes/public/SectionValues.m (revision 26) +++ /issm/trunk/src/m/classes/public/SectionValues.m (revision 27) @@ -83,6 +83,6 @@ %Get bed and surface for each 2d point, offset to make sure that it is inside the glacier system offset=10^-10; - bed=griddata_mesh_to_mesh_3d(md.elements,md.x,md.y,md.z,md.bed,X,Y,Z,'node')+offset; - surface=griddata_mesh_to_mesh_3d(md.elements,md.x,md.y,md.z,md.surface,X,Y,Z,'node')-offset; + bed=griddata_mesh_to_mesh(md.elements2d,md.x2d,md.y2d,project2d(md,md.bed,1),X,Y,'node')+offset; + surface=griddata_mesh_to_mesh(md.elements2d,md.x2d,md.y2d,project2d(md,md.surface,1),X,Y,'node')+offset; %Some useful parameters Index: /issm/trunk/src/m/classes/public/ThicknessCorrection.m =================================================================== --- /issm/trunk/src/m/classes/public/ThicknessCorrection.m (revision 26) +++ /issm/trunk/src/m/classes/public/ThicknessCorrection.m (revision 27) @@ -20,6 +20,6 @@ in=ArgusContourToMesh(md.elements,md.x,md.y,expread(filename,1),'node',1); -pos_shelf=find(in & md.gridoniceshelf); -pos_sheet=find(in & ~md.gridoniceshelf); +pos_shelf=find(in & ~md.gridonicesheet); +pos_sheet=find(in & md.gridonicesheet); for i=1:length(pos_shelf) %search the grid on ice sheet the closest to i Index: /issm/trunk/src/m/classes/public/displaycontrol.m =================================================================== --- /issm/trunk/src/m/classes/public/displaycontrol.m (revision 26) +++ /issm/trunk/src/m/classes/public/displaycontrol.m (revision 27) @@ -10,5 +10,16 @@ disp(sprintf(' ''%s''','control')); -disp(sprintf(' control_type: %s (control type, ex: ''drag'', or ''B'')',md.control_type)); +%control type +control_string=''; +for i=1:length(md.control_type), + parameter=md.control_type{i}; + %check this parameter is a field from model! + if ~isfield(struct(md),parameter), + error('displaysolutionparameters error message: one of the control type parameters does not exist!'); + end + control_string=[control_string parameter ' and ']; +end +control_string=control_string(1:length(control_string)-5); +disp(sprintf(' control_type: %s %s',control_string,'(list of parameters where inverse control is carried out; ex: {''drag''}, or {''drag'',''B''})')); disp(sprintf(' fit: (%i) (''absolute: 0'', ''relative: 1'', or ''logarithmic: 2''. default is ''absolute: 0'', for each optimization steps)',length(md.fit))); disp(sprintf(' meanvel: %g (velocity scaling factor when evaluating relative or logarithmic misfit)',md.meanvel)); Index: /issm/trunk/src/m/classes/public/displayqmu.m =================================================================== --- /issm/trunk/src/m/classes/public/displayqmu.m (revision 26) +++ /issm/trunk/src/m/classes/public/displayqmu.m (revision 27) @@ -10,4 +10,5 @@ disp(sprintf(' ''%s''','qmu using Dakota')); + disp(sprintf(' variables: (arrays of each variable class)')); fnames=fieldnames(md.variables); @@ -16,4 +17,5 @@ fnames{i},size(md.variables.(fnames{i})),class(md.variables.(fnames{i})))); end + disp(sprintf(' responses: (arrays of each response class)')); fnames=fieldnames(md.responses); @@ -22,16 +24,32 @@ fnames{i},size(md.responses.(fnames{i})),class(md.responses.(fnames{i})))); end + +disp(sprintf(' qmu_method: (array of dakota_method class)')); +for i=1:numel(md.qmu_method); + if strcmp(class(md.qmu_method(i)),'dakota_method') + disp(sprintf(' method%s : ''%s''',... + string_dim(md.qmu_method,i),md.qmu_method(i).method)); + end +end + +for i=1:numel(md.qmu_params) + disp(sprintf(' qmu_params%s: (array of method-independent parameters)',... + string_dim(md.qmu_params,i))); + fnames=fieldnames(md.qmu_params(i)); + maxlen=0; + for j=1:numel(fnames) + maxlen=max(maxlen,length(fnames{j})); + end + + for j=1:numel(fnames) + disp(sprintf([' %-' num2str(maxlen+1) 's: %s'],... + fnames{j},any2str(md.qmu_params(i).(fnames{j})))); + end +end + disp(sprintf(' dakotaresults: 1x%i {dvar,rfunc,scm,pcm,srcm,prcm}',length(md.dakotaresults))); if isempty(md.dakotain), disp(sprintf(' dakotain: N/A')); else disp(sprintf(' dakotain: not displayed (can be accessed by typing md.dakotain)'));end if isempty(md.dakotaout), disp(sprintf(' dakotaout: N/A')); else disp(sprintf(' dakotaout: not displayed (can be accessed by typing md.dakotaout)'));end if isempty(md.dakotadat), disp(sprintf(' dakotadat: N/A')); else disp(sprintf(' dakotadat: not displayed (can be accessed by typing md.dakotadat)'));end -disp(sprintf(' method: ''%s'' (''local reliability'' or ''sampling'')',md.method)); -disp(sprintf(' seed: %g (random seed number)',md.seed)); -disp(sprintf(' samples: %i (number of samples)',md.samples)); -disp(sprintf(' sample_type: ''%s'' (type of sampling, ''random'' or ''lhs'')',md.sample_type)); -disp(sprintf(' interval_type: ''%s'' (FD ''forward'' or ''central'')',md.interval_type)); -disp(sprintf(' fd_gradient_step_size: %g (FD gradient step size, default= .001) ',md.fd_gradient_step_size)); -disp(sprintf(' evaluation_concurrency: %i (evaluation concurrency level)',md.evaluation_concurrency)); -disp(sprintf(' npart: %i (number of partitions for semi-descrete qmu)',md.npart)); -disp(sprintf(' analysis_driver : ''%s'' (''matlab'' or name of driver for Dakota analysis)',md.analysis_driver)); -disp(sprintf(' analysis_components: ''%s'' (Matlab script file for Matlab driver)',md.analysis_components)); +disp(sprintf(' npart : %i (number of partitions for semi-descrete qmu)',md.npart)); +disp(sprintf(' numrifts: %i (number of rifts for semi-descrete qmu)',md.numrifts)); Index: /issm/trunk/src/m/classes/public/extrude.m =================================================================== --- /issm/trunk/src/m/classes/public/extrude.m (revision 26) +++ /issm/trunk/src/m/classes/public/extrude.m (revision 27) @@ -152,5 +152,5 @@ md.segmentonneumann_diag=segmentonneumann_diag; md.gridondirichlet_prog=project3d(md,md.gridondirichlet_prog,'node'); -md.dirichletvalues_prog=project3d(md,md.dirichletvalues_prog,'node',1); +md.dirichletvalues_prog=project3d(md,md.dirichletvalues_prog,'node'); %md.segmentonneumann_prog=[tproj(md.segmentonneumann_prog(:,1)) tproj(md.segmentonneumann_prog(:,2)) tproj2d_el(md.segmentonneumann_prog(:,3))]; %md.segmentonneumann_prog2=[tproj(md.segmentonneumann_prog2(:,1)) tproj(md.segmentonneumann_prog2(:,2)) tproj2d_el(md.segmentonneumann_prog2(:,3))]; Index: /issm/trunk/src/m/classes/public/ismodelselfconsistent.m =================================================================== --- /issm/trunk/src/m/classes/public/ismodelselfconsistent.m (revision 26) +++ /issm/trunk/src/m/classes/public/ismodelselfconsistent.m (revision 27) @@ -5,18 +5,313 @@ % bool=ismodelselfconsistent(md,solutiontype,package), -if strcmpi(solutiontype,'mesh'), - %this solution is a little special. It should come right after the md=model; operation. So a lot less checks! - - bool=1; - return; -end - %tolerance we use in litmus tests for the consistency of the model tolerance=10^-12; if (nargin~=3 ) - IsModelSelfConsistentUsage(); - error(' '); -end - + help ismodelselfconsistent + error('ismodelselfconsistent error message: wrong usage'); +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% COMMON CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%COUNTER +if md.counter<3, + disp(['model ' md.name ' is not correctly configured. You forgot one step in the following sequence (mesh, geography, parameterize)!']); + bool=0;return; +end + +%NAME +if isempty(md.name), + disp(['model is not correctly configured: missing name!']); + bool=0;return; +end + +%MESH +if md.numberofelements<=0, + disp(['model ' md.name ' does not have any elements!']); + bool=0; return; +end +if md.numberofgrids<=0, + disp(['model ' md.name ' does not have any grids!']); + bool=0; return; +end + +%ELEMENTSTYPE +if size(md.elements_type,1)~=md.numberofelements | size(md.elements_type,2)~=2, + disp(['Types of elements have not been set properly, run setelementstype first']) + bool=0;return; +end +if any(ones(md.numberofelements,1)-((md.elements_type(:,1)==hutterenum) + (md.elements_type(:,1)==macayealenum) + (md.elements_type(:,1)==pattynenum))) + disp(['Types of elements have not been set properly, run setelementstype first']) + bool=0;return; +end +if any(ones(md.numberofelements,1)-((md.elements_type(:,2)==stokesenum) + (md.elements_type(:,2)==noneenum))) + disp(['Types of elements have not been set properly, run setelementstype first']) + bool=0;return; +end +if strcmpi(md.type,'2d'), + if (ismember(pattynenum,md.elements_type(:,1)) | ismember(stokesenum,md.elements_type(:,2))), + disp(['For a 2d model, only MacAyeal''s and Hutter''s elements are allowed']); + bool=0;return; + end +end + +%NO NAN +fields={'numberofelements','numberofgrids','x','y','z','segmentonneumann_diag','drag','drag_type','p','q','gridondirichlet_diag',... + 'rho_ice','rho_water','B','elementoniceshelf','surface','thickness','bed','g','lowmem','sparsity','nsteps','maxiter',... + 'tolx','np','eps_rel','exclusive','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'}; +for i=1:length(fields), + if ~isempty(eval(['md.' char(fields(i))])), + if find(isnan(eval(['md.' char(fields(i))]))), + disp(['model ' md.name ' has an NaN value in field ' char(fields(i)) '!']); + bool=0; return; + end + end +end + +%FIELDS > 0 +fields={'numberofelements','numberofgrids','elements','segmentonneumann_diag','drag','drag_type','p','q','gridondirichlet_diag',... + 'rho_ice','rho_water','B','elementoniceshelf','thickness','g','eps_rel','eps_abs','nsteps','maxiter','tolx','exclusive',... + 'sparsity','lowmem','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'}; +for i=1:length(fields), + if ~isempty(eval(['md.' char(fields(i))])), + if find((eval(['md.' char(fields(i))]))<0), + disp(['model ' md.name ' has a <0 value in field ' char(fields(i)) '!']); + bool=0; return; + end + end +end +if any(md.p<=0), + disp(['model ' md.name ' has some p<0 friction coefficientsin sliding law']); + bool=0; return; +end + +%FIELDS ~=0 +fields={'numberofelements','numberofgrids','elements','segmentonneumann_diag','drag_type',... + 'rho_ice','rho_water','B','thickness','g','eps_rel','eps_abs','maxiter','tolx',... + 'sparsity','deltaH','DeltaH','timeacc','timedec'}; +for i=1:length(fields), + if ~isempty(eval(['md.' char(fields(i))])), + if find((eval(['md.' char(fields(i))]))==0), + disp(['model ' md.name ' has a =0 value in field ' char(fields(i)) '!']); + bool=0; return; + end + end +end + +%SIZE NUMBEROFELEMENTS +fields={'elements','p','q','elementoniceshelf','n','elementonbed'}; +for i=1:size(fields,2), + if (size(eval(['md.' char(fields(i))]),1)~=md.numberofelements), + disp(['model ' md.name ' field ' char(fields(i)) ' should be of size ' num2str(md.numberofelements) '!']); + bool=0; return; + end +end + +%SIZE NUMBEROFGRIDS +fields={'x','y','z','B','drag','gridondirichlet_diag','melting','accumulation','surface','thickness','bed','gridonbed','gridonsurface'}; +for i=1:length(fields), + if find(length(eval(['md.' char(fields(i))]))~=md.numberofgrids), + disp(['model ' md.name ' field ' char(fields(i)) ' should be of size ' num2str(md.numberofgrids) '!']); + bool=0; return; + end +end + +%THICKNESS = SURFACE - BED +if find((md.thickness-md.surface+md.bed)>tolerance), + disp(['model ' md.name ' violates the equality thickness=surface-bed!']); + bool=0; return; +end + +%RIFTS +if md.numrifts, + if ~strcmpi(md.type,'2d'), + disp(['Models with rifts are only supported in 2d for now!']); + bool=0;return; + end +end +if ~isstruct(md.rifts), + if ~isempty(find(md.segmentmarkers>=2)), + %We have segments with rift markers, but no rift structure! + disp(['model ' md.name ' should be processed for rifts (run meshprocessrifts)!']); + bool=0; return; + end +end + +%ARTIFICIAL DIFFUSIVITY +if ~isscalar(md.artificial_diffusivity), + disp('artificial_diffusivity should be a scalar (1 or 0)'); + bool=0;return; +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SOLUTION CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%DIAGNOSTIC +if strcmpi(solutiontype,'diagnostic') + + %HUTTER ON ICESHELF WARNING + if any(md.elements_type(:,1)==hutterenum & md.elementoniceshelf), + disp(sprintf('\n !!! Warning: Hutter''s model is not consistent on ice shelves !!!\n')); + end + + %SINGULAR + if ~any(md.gridondirichlet_diag), + disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!']) + bool=0;return; + end + + %DIRICHLET VALUES + if isempty(md.gridondirichlet_diag) | isempty(md.dirichletvalues_diag), + disp(['model ' md.name ' is not well posed. Missing dirichlet values for diagnostic run']); + bool=0;return; + end + + %DIRICHLET IF THICKNESS <= 0 + if ~isempty(find(md.thickness<=0)), + pos=find(md.thickness<=0); + if ~isempty(find(md.gridondirichlet_diag(pos)==0)), + disp(['model ' md.name ' has some grids with 0 thickness']); + bool=0; return; + end + end +end + +%PROGNOSTIC +if strcmp(solutiontype,'prognostic'), + %old testing + % if isempty(find(md.gridondirichlet_diag)), + % disp(['model ' md.name ' diagnostic is not well posed (singular). You need at least one grid with fixed velocity!']) + % bool=0;return; + % end + % if isempty(find(md.gridondirichlet_prog)), + % disp(['model ' md.name ' prognostic is not well posed (singular). You need at least one grid with fixed thickness!']) + % bool=0;return; + % end + + %VELOCITIES + if (size(md.vx,1)~=md.numberofgrids | size(md.vy,1)~=md.numberofgrids), + disp(['a 3d velocity is required. Run ''diagnostic'' solution first!']) + bool=0; return; + end +end + +%THERMAL TRANSIENT +if strcmp(solutiontype,'thermaltransient') + + %INITIAL TEMPERATURE, MELTING AND ACCUMULATION + if isempty(md.temperature), + disp(['An initial temperature is needed for a transient thermal computation']) + bool=0;return; + end + if isstruct(md.temperature) | isstruct(md.melting) | isstruct(md.accumulation), + disp(['The initial temperature, melting or accumulation should be a list and not a structure']) + bool=0;return; + end +end + +%THERMAL STEADY AND THERMAL TRANSIENT +if strcmpi(solutiontype,'thermalsteady') | strcmp(solutiontype,'thermaltransient'), + + %EXTRUSION + if strcmp(md.type,'2d'), + disp(['For a ' solutiontype ' computation, the model must be 3d, extrude it first!']) + bool=0;return; + end + + %VELOCITIES AND PRESSURE + if (isempty(md.vx) | isnan(md.vx) | isempty(md.vy) | isnan(md.vy) | isempty(md.vz) | isnan(md.vz)), + disp(['a 3d velocity is required. Run ''diagnostic'' solution first!']) + bool=0;return; + end + if (isempty(md.pressure) | isnan(md.pressure)), + disp(['pressure is required. Run ''diagnostic'' solution first!']) + bool=0;return; + end +end + +%THERMAL TRANSIENT AND TRANSIENT +if strcmp(solutiontype,'thermaltransient') | strcmp(solutiontype,'transient'), + + %DT and NDT + fields={'dt','ndt'}; + for i=1:length(fields), + if find((eval(['md.' char(fields(i))]))<0), + disp(['model ' md.name ' has a <0 value in field ' char(fields(i)) '!']); + bool=0; return; + end + end + + %INITIAL TEMPERATURE + if isstruct(md.temperature), + disp(['The initial temperature should be empty or a list but not a structure']) + bool=0;return; + end +end + +%PARAMETERS +if strcmp(solutiontype,'parameters') + + %OUTPUT + if ~iscell(md.parameteroutput) + disp(['parameteroutput field must be a cell, example {''strainrate'',''stress'',''deviatoricstress'',''viscousheating''}']); + bool=0; return; + end + for i=1:length(md.parameteroutput) + if ~strcmpi(md.parameteroutput(i),'strainrate') & ~strcmpi(md.parameteroutput(i),'stress') & ~strcmpi(md.parameteroutput(i),'deviatoricstress') & ~strcmpi(md.parameteroutput(i),'viscousheating') ... + & ~strcmpi(md.parameteroutput(i),'pressure_elem') & ~strcmpi(md.parameteroutput(i),'stress_bed') & ~strcmpi(md.parameteroutput(i),'stress_surface') + disp(['one of the parameteroutput is not supported yet']); + bool=0; return; + end + end + %VELOCITY + if ~(size(md.vx,1)==md.numberofgrids & size(md.vy,1)==md.numberofgrids & (size(md.vz,1)==md.numberofgrids | strcmpi(md.type,'2d'))) + disp(['velocities are required!']); + bool=0;return; + end + + %ACCELERATION + if strcmpi(md.type,'2d') & md.acceleration + disp(['solution requires acceleration=0']); + bool=0;return; + end + + %HUTTER + if any(md.elements_type(:,1)==hutterenum); + disp(['The model has Hutter''s elements. Impossible to compute parameters']); + bool=0;return; + end +end + +%CONTROL +if strcmpi(solutiontype,'control'), + + %CONTROL TYPE + if ~iscell(md.control_type), + disp('control_type should be a cell array of strings'); + bool=0;return; + end + + %LENGTH CONTROL FIELDS + if (length(md.maxiter)~=md.nsteps | length(md.optscal)~=md.nsteps | length(md.fit)~=md.nsteps) + disp('maxiter, optscal and fit must have the length specified by nsteps') + bool=0;return; + end + + %FIT + if sum((double(md.fit==1) + double(md.fit==0) + double(md.fit==2))==1)~=md.nsteps + disp('wrong fits: fit should be a vector of size nsteps holding 0, 1 and 2 only') + bool=0;return; + end + + %OBSERVED VELOCITIES + fields={'vx_obs','vy_obs'}; + for i=1:length(fields), + if find(length(eval(['md.' char(fields(i))]))~=md.numberofgrids), + disp(['model ' md.name ' field ' char(fields(i)) ' should be of size ' num2str(md.numberofgrids) '!']); + bool=0; return; + end + end +end + +%QMU if strcmpi(solutiontype,'qmu'), if ~strcmpi(md.cluster,'none'), @@ -28,341 +323,95 @@ end -if isempty(md.name), - disp(['model is not correctly configured: missing name!']); - bool=0;return; -end - -if md.counter<3, - disp(['model ' md.name ' is not correctly configured. You forgot one step in the following sequence (mesh, geography, parameterize)!']); - bool=0;return; -end - -if md.numberofelements<=0, - disp(['model ' md.name ' does not have any elements!']); - bool=0; return; -end - -if md.numberofgrids<=0, - disp(['model ' md.name ' does not have any grids!']); - bool=0; return; -end - -%check elementstype -if size(md.elements_type,1)~=md.numberofelements | size(md.elements_type,2)~=2, - disp(['Types of elements have not been set properly, run setelementstype first']) - bool=0;return; -end -if any(ones(md.numberofelements,1)-((md.elements_type(:,1)==hutterenum) + (md.elements_type(:,1)==macayealenum) + (md.elements_type(:,1)==pattynenum))) - disp(['Types of elements have not been set properly, run setelementstype first']) - bool=0;return; -end -if any(ones(md.numberofelements,1)-((md.elements_type(:,2)==stokesenum) + (md.elements_type(:,2)==noneenum))) - disp(['Types of elements have not been set properly, run setelementstype first']) - bool=0;return; -end -if strcmpi(md.type,'2d'), - if (ismember(pattynenum,md.elements_type(:,1)) | ismember(stokesenum,md.elements_type(:,2))), - disp(['For a 2d model, only MacAyeal''s and Hutter''s elements are allowed']); - bool=0;return; - end -end -%check that if there are rifts, model is 2d -if md.numrifts, - if ~strcmpi(md.type,'2d'), - disp(['Models with rifts are only supported in 2d for now!']); - bool=0;return; - end -end - -%If there are rifts, check that the solver is set to mumps in parallel -if md.numrifts, - if ~strcmpi(md.cluster,'none') - if isempty(findstr('aijmumps',md.solverstring)), - disp(['For a 2d model with rifts, running on a cluster, you should use a direct solver like MUMPS!']); +%MESH +if strcmpi(solutiontype,'mesh'), + %this solution is a little special. It should come right after the md=model; operation. So a lot less checks! + + bool=1; + return; +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PACKAGE CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%CIELO +if strcmpi(package,'cielo'), + + %NAN VALUES + fields={'time','sparsity'}; + for i=1:length(fields), + if ~isempty(eval(['md.' char(fields(i))])), + if find(isnan(eval(['md.' char(fields(i))]))), + disp(['model ' md.name ' has an NaN value in field ' char(fields(i)) '!']); + bool=0; return; + end + end + end + + %FIELD > 0 + fields={'time','sparsity'}; + for i=1:length(fields), + if ~isempty(eval(['md.' char(fields(i))])), + if find((eval(['md.' char(fields(i))]))<0), + disp(['model ' md.name ' has a <0 value in field ' char(fields(i)) '!']); + bool=0; return; + end + end + end + + %FIELD ~= 0 + fields={'time','sparsity'}; + for i=1:length(fields), + if ~isempty(eval(['md.' char(fields(i))])), + if find((eval(['md.' char(fields(i))]))==0), + disp(['model ' md.name ' has a =0 value in field ' char(fields(i)) '!']); + bool=0; return; + end + end + end + + %SPARSITY BETWEEN 0 AND 1 + if ( (md.sparsity<=0) | (md.sparsity>1)), + disp(['model ' md.name ' sparsity should be inside the ]0 1] range']); + bool=0; return; + end + + %RIFTS + if md.numrifts, + if ~strcmpi(md.cluster,'none') + if isempty(findstr('aijmumps',md.solverstring)), + disp(['For a 2d model with rifts, running on a cluster, you should use a direct solver like MUMPS!']); + bool=0;return; + end + end + end + + %CONNECTIVITY + if strcmpi(md.type,'2d'), + if md.connectivity<9, + disp('connectivity should be at least 9 for 2d models'); bool=0;return; end end -end - -%check elementstype in diagnostic -if strcmpi(solutiontype,'diagnostic') - if any(md.elements_type(:,1)==hutterenum & md.elementoniceshelf), - disp(sprintf('\n !!! Warning: Hutter''s model is not consistent on ice shelves !!!\n')); - end -end - -%Check that mesh is 3d -if strcmpi(solutiontype,'thermalsteady') | strcmp(solutiontype,'thermaltransient'), - if strcmp(md.type,'2d'), - disp(['For a ' solutiontype ' computation, the model must be 3d, extrude it first!']) - bool=0;return; - end -end - -%check artificial_diffusivity is an integer -if ~isscalar(md.artificial_diffusivity), - disp('artificial_diffusivity should be a scalar (1 or 0)'); - bool=0;return; -end - -%check on connectivity -if strcmpi(md.type,'2d'), - if md.connectivity<9, - disp('connectivity should be at least 9 for 2d models'); - bool=0;return; - end -end - -if strcmpi(md.type,'3d'), - if md.connectivity<24, - disp('connectivity should be at least 24 for 3d models'); - bool=0;return; - end -end - - - - - -%Check previous computation -if strcmp(solutiontype,'thermaltransient') | strcmp(solutiontype,'thermalsteady') - if (isempty(md.vx) | isempty(md.vy) | isempty(md.vz)), - disp(['a 3d velocity is required. Run ''diagnostic'' solution first!']) - bool=0;return; - end -end - -if strcmp(solutiontype,'thermaltransient') | strcmp(solutiontype,'thermalsteady') - if (isempty(md.pressure) | isnan(md.pressure)), - disp(['pressure is required. Run ''diagnostic'' solution first!']) - bool=0;return; - end - - if (isempty(md.vx) | isnan(md.vx) | isempty(md.vy) | isnan(md.vy) | isempty(md.vz) | isnan(md.vz)), - disp(['velocities are required!']); - bool=0;return; - end -end - -if strcmp(solutiontype,'thermaltransient') - if isempty(md.temperature), - disp(['An initial temperature is needed for a transient thermal computation']) - bool=0;return; - end - if isstruct(md.temperature) | isstruct(md.melting), - disp(['The initial temperature or melting should be a list and not a structure']) - bool=0;return; - end - -end - -if strcmp(solutiontype,'transient') - if isstruct(md.temperature), - disp(['The initial temperature should be empty or a list but not a structure']) - bool=0;return; - end - -end - -if strcmp(solutiontype,'parameters') - if ~(size(md.vx,1)==md.numberofgrids & size(md.vy,1)==md.numberofgrids & (size(md.vz,1)==md.numberofgrids | strcmpi(md.type,'2d'))) - disp(['velocities are required!']); - bool=0;return; - end - if strcmpi(md.type,'2d') & md.acceleration - disp(['solution requires acceleration=0']); - bool=0;return; - end - if any(md.elements_type(:,1)==hutterenum); - disp(['The model has Hutter''s elements. Impossible to compute parameters']); - bool=0;return; - end -end - - -%Check that problem is not singular -if strcmp(solutiontype,'diagnostic'), - if isempty(find(md.gridondirichlet_diag)), - disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!']) - bool=0;return; - end -end - -%Check we have values for dirichletvalues_diag -if strcmp(solutiontype,'diagnostic'), - if isempty(md.gridondirichlet_diag) | isempty(md.dirichletvalues_diag), - disp(['model ' md.name ' is not well posed. Missing dirichlet values for diagnostic run']); - bool=0;return; - end -end - - - -if strcmp(solutiontype,'prognostic'), - %old testing -% if isempty(find(md.gridondirichlet_diag)), -% disp(['model ' md.name ' diagnostic is not well posed (singular). You need at least one grid with fixed velocity!']) -% bool=0;return; -% end -% if isempty(find(md.gridondirichlet_prog)), -% disp(['model ' md.name ' prognostic is not well posed (singular). You need at least one grid with fixed thickness!']) -% bool=0;return; -% end - if (size(md.vx,1)~=md.numberofgrids | size(md.vy,1)~=md.numberofgrids), - disp(['a 3d velocity is required. Run ''diagnostic'' solution first!']) - bool=0; return; - end -end -%some checks for control methods -if strcmp(solutiontype,'control') - %check fields used by control - if (length(md.maxiter)~=md.nsteps | length(md.optscal)~=md.nsteps | length(md.fit)~=md.nsteps) - disp('maxiter, optscal and fit must have the length specified by nsteps') - bool=0;return; - end - if sum((double(md.fit==1) + double(md.fit==0) + double(md.fit==2))==1)~=md.nsteps - disp('wrong fits: fit should be a vector of size nsteps holding 0, 1 and 2 only') - bool=0;return; - end -end - - -%Check that no NaN values popped up: -fields={'numberofelements','numberofgrids','x','y','z','segmentonneumann_diag','drag','drag_type','p','q','gridondirichlet_diag',... - 'rho_ice','rho_water','B','elementoniceshelf','surface','thickness','bed','g','lowmem','sparsity','nsteps','maxiter',... - 'tolx','np','time','eps_rel','exclusive','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'}; -for i=1:length(fields), - if ~isempty(eval(['md.' char(fields(i))])), - if find(isnan(eval(['md.' char(fields(i))]))), - disp(['model ' md.name ' has an NaN value in field ' char(fields(i)) '!']); - bool=0; return; - end - end -end - -%Check that these fields are > 0 -fields={'numberofelements','numberofgrids','elements','segmentonneumann_diag','drag','drag_type','p','q','gridondirichlet_diag',... - 'rho_ice','rho_water','B','elementoniceshelf','thickness','g','eps_rel','eps_abs','nsteps','maxiter','tolx','np','time','exclusive',... - 'sparsity','lowmem','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'}; -for i=1:length(fields), - if ~isempty(eval(['md.' char(fields(i))])), - if find((eval(['md.' char(fields(i))]))<0), - disp(['model ' md.name ' has a <0 value in field ' char(fields(i)) '!']); - bool=0; return; - end - end -end - -%Check that these fields are ~=0 -fields={'numberofelements','numberofgrids','elements','segmentonneumann_diag','drag_type',... - 'rho_ice','rho_water','B','thickness','g','eps_rel','eps_abs','maxiter','tolx','np','time',... - 'sparsity','deltaH','DeltaH','timeacc','timedec'}; -for i=1:length(fields), - if ~isempty(eval(['md.' char(fields(i))])), - if find((eval(['md.' char(fields(i))]))==0), - disp(['model ' md.name ' has a =0 value in field ' char(fields(i)) '!']); - bool=0; return; - end - end -end - -%Check that these fields are of size md.numberofelements -fields={'elements','p','q','elementoniceshelf','n','elementonbed'}; -for i=1:size(fields,2), - if (size(eval(['md.' char(fields(i))]),1)~=md.numberofelements), - disp(['model ' md.name ' field ' char(fields(i)) ' should be of size ' num2str(md.numberofelements) '!']); - bool=0; return; - end -end - -%Check that these fields are of size md.numberofgrids -fields={'melting','x','y','z','B','drag','gridondirichlet_diag','vx_obs','vy_obs','surface','thickness','bed','gridonbed','gridonsurface'}; -for i=1:length(fields), - if find(length(eval(['md.' char(fields(i))]))~=md.numberofgrids), - disp(['model ' md.name ' field ' char(fields(i)) ' should be of size ' num2str(md.numberofgrids) '!']); - bool=0; return; - end -end - -%Make sure lowmem is either 0 or 1 -if ((md.lowmem ~= 1) & (md.lowmem~=0)), - disp(['model ' md.name ' lowmem field should be 0 or 1']); - bool=0; return; -end - -%Make sure sparsity is between 0 and 1: -if ( (md.sparsity<=0) | (md.sparsity>1)), - disp(['model ' md.name ' sparsity should be inside the ]0 1] range']); - bool=0; return; -end - -%Make sure thickness=surface-bed; -if find((md.thickness-md.surface+md.bed)>tolerance), - disp(['model ' md.name ' violates the equality thickness=surface-bed!']); - bool=0; return; -end - -%Check that rifts were processed correctly: -if ~isstruct(md.rifts), - if ~isempty(find(md.segmentmarkers>=2)), - %We have segments with rift markers, but no rift structure! - disp(['model ' md.name ' should be processed for rifts (run meshprocessrifts)!']); - bool=0; return; - end -end - - -%Check that thickness is not null. If is it, then the corresponding grids should be spc'd. -if ~isempty(find(md.thickness<=0)), - pos=find(md.thickness<=0); - if ~isempty(find(md.gridondirichlet_diag(pos)==0)), - disp(['model ' md.name ' has some grids with 0 thickness']); - bool=0; return; - end -end - -%Check that p>0 in the sliding law u = k * sigma ^ p * Neff ^ q : -if ~isempty(find(md.p<=0)), - disp(['model ' md.name ' has some p<0 friction coefficientsin sliding law']); - bool=0; return; -end - - -%Check parameteroutput -if ~iscell(md.parameteroutput) - disp(['parameteroutput field must be a cell, example {''strainrate'',''stress'',''deviatoricstress'',''viscousheating''}']); - bool=0; return; -end - - -for i=1:length(md.parameteroutput) - if ~strcmpi(md.parameteroutput(i),'strainrate') & ~strcmpi(md.parameteroutput(i),'stress') & ~strcmpi(md.parameteroutput(i),'deviatoricstress') & ~strcmpi(md.parameteroutput(i),'viscousheating') ... - & ~strcmpi(md.parameteroutput(i),'pressure_elem') & ~strcmpi(md.parameteroutput(i),'stress_bed') & ~strcmpi(md.parameteroutput(i),'stress_surface') - disp(['one of the parameteroutput is not supported yet']); - bool=0; return; - end -end - -%check parallel settings -if ~strcmpi(md.cluster,'none'), - if ~strcmpi(package,'cielo'), - disp(['parallel runs only allowed for ''cielo'' package']); - bool=0; return; - end -end - -%Put as many checks as you want here: - + if strcmpi(md.type,'3d'), + if md.connectivity<24, + disp('connectivity should be at least 24 for 3d models'); + bool=0;return; + end + end + + %NP + if md.np==0, + disp(['model ' md.name ' has a =0 value in field ' md.np ' !']); + elseif md.np<0, + disp(['model ' md.name ' has a negative value in field ' md.np ' !']); + end + + %LOWMEM = 0 or 1 + if ((md.lowmem ~= 1) & (md.lowmem~=0)), + disp(['model ' md.name ' lowmem field should be 0 or 1']); + bool=0; return; + end +end %No problems, just return 1; bool=1; return; - -function IsModelSelfConsistentUsage() -disp(' '); -disp(' IsModelSelfConsistent usage: flag=IsModelSelfConsistent(model)'); -disp(' where model is an instance of the @model class, and flag a boolean'); - - - - Index: /issm/trunk/src/m/classes/public/isresultconsistent.m =================================================================== --- /issm/trunk/src/m/classes/public/isresultconsistent.m (revision 26) +++ /issm/trunk/src/m/classes/public/isresultconsistent.m (revision 27) @@ -21,5 +21,5 @@ if ~isnan(md.penalties), for i=1:size(md.penalties,1) - for n=2:md.numlayers + for n=2:size(md.penalties,2) relativex=(md.vx(md.penalties(i,1))-md.vx(md.penalties(i,n)))./(md.vx(md.penalties(i,1))); if md.vx(md.penalties(i,1))==md.vx(md.penalties(i,n)), relativex=0; end @@ -74,5 +74,5 @@ if ~isnan(md.penalties), for i=1:size(md.penalties,1) - for n=2:md.numlayers + for n=2:size(md.penalties,2) relativex=(md.transient_results(iter).vx(md.penalties(i,1))-md.transient_results(iter).vx(md.penalties(i,n)))./(md.transient_results(iter).vx(md.penalties(i,1))); if md.transient_results(iter).vx(md.penalties(i,1))==md.transient_results(iter).vx(md.penalties(i,n)), relativex=0; end Index: sm/trunk/src/m/classes/public/mesh2.m =================================================================== --- /issm/trunk/src/m/classes/public/mesh2.m (revision 26) +++ (revision ) @@ -1,73 +1,0 @@ -function md2=mesh2(md,domainname,varargin) -%MESH - create model mesh -% -% This routine creates a model mesh using TriMesh and a domain outline, to within a certain resolution -% where md is a @model object, domainname is the name of an Argus domain outline file, -% and resolution is a characteristic length for the mesh (same unit as the domain outline -% unit). Riftname is an optional argument (Argus domain outline) describing rifts. -% -% Usage: -% md2=mesh2(md,domainname,resolution) -% or md2=mesh2(md,domainname,riftname, resolution) -% -% Examples: -% md2=mesh2(md,'DOmainOutline.exp',1000); -% md2=mesh2(md,'DOmainOutline.exp','Rifts.exp',1500); - -%Figure out a characteristic area. Resolution is a grid oriented concept (ex a 1000m resolution grid would -%be made of 1000*1000 area squares). -if (nargin==3), - resolution=varargin{1}; - riftname=''; -end -if (nargin==4), - riftname=varargin{1}; - resolution=varargin{2}; -end - -%Check that mesh was not already run, and warn user: -if subsref(md,struct('type','.','subs','counter'))>=1, - choice=input('This model already has a mesh. Are you sure you want to go ahead? (y/n)','s'); - if ~strcmp(choice,'y') - error('no meshing done ... exiting'); - end -end - -area=resolution^2; -%Initialize return model; -md2=model; - -%Mesh using TriMesh -if strcmp(riftname,''), - [elements,x,y,segments,segmentmarkers]=TriMesh(domainname,area,'yes'); -else - [elements,x,y,segments,segmentmarkers]=TriMesh(domainname,riftname,area,'yes'); -end - -md2=subsasgn(md2,struct('type','.','subs','elements'),elements); -md2=subsasgn(md2,struct('type','.','subs','segmentmarkers'),segmentmarkers); -md2=subsasgn(md2,struct('type','.','subs','segments'),segments); -md2=subsasgn(md2,struct('type','.','subs','y'),y); -md2=subsasgn(md2,struct('type','.','subs','x'),x); - -%Fill in rest of fields: -md2=subsasgn(md2,struct('type','.','subs','numberofelements'),length(elements)); -md2=subsasgn(md2,struct('type','.','subs','numberofgrids'),length(x)); -md2=subsasgn(md2,struct('type','.','subs','z'),zeros(length(x),1)); -gridonboundary=zeros(length(x),1); gridonboundary(segments(:,1:2))=1; -md2=subsasgn(md2,struct('type','.','subs','gridonboundary'),gridonboundary); -md2=subsasgn(md2,struct('type','.','subs','elements_type'),3*ones(md2.numberofelements,1)); % type determined by number of grids per element - -%outline names -md2=subsasgn(md2,struct('type','.','subs','domainoutline'),readfile(domainname)); -if strcmp(riftname,''), - md2=subsasgn(md2,struct('type','.','subs','riftoutline'),''); -else - md2=subsasgn(md2,struct('type','.','subs','riftoutline'),readfile(riftname)); -end - -%type of model -md2=subsasgn(md2,struct('type','.','subs','type'),'2d'); - -%augment counter keeping track of what has been done to this model -md2=subsasgn(md2,struct('type','.','subs','counter'),1); Index: /issm/trunk/src/m/classes/public/parameterize.m =================================================================== --- /issm/trunk/src/m/classes/public/parameterize.m (revision 26) +++ /issm/trunk/src/m/classes/public/parameterize.m (revision 27) @@ -27,21 +27,22 @@ end -%Try and run parameter file. We try to capture the error message, but this is not allowed for -%lower versions of matlab. +%Try and run parameter file. +eval(['system( '' cp ' parametername ' TemporaryParameterFile.m'');']); -stringfile=readfile(parametername); -stringlines=strsplit(stringfile,char(10)); %split into lines - -for i=1:length(stringlines), - a=stringlines{i}; - try, - eval(stringlines{i}); - catch me, - disp(' '); - disp(['Error in ==> ' parametername ' at ' num2str(i)]); - disp(me.message); - disp(stringlines{i}); - error(' '); +try, + TemporaryParameterFile + system('rm TemporaryParameterFile.m'); +catch me, + system('rm TemporaryParameterFile.m'); + me2=struct('message',me.message,'stack',me.stack); + for i=1:length(me2.stack) + if (length(me2.stack(i).file)>=24 & strcmp(me2.stack(i).file(end-23:end),'TemporaryParameterFile.m')) + me2.stack(i).file=[me2.stack(i).file(1:end-24) parametername]; + end + if strcmp(me2.stack(i).name,'TemporaryParameterFile'), + me2.stack(i).name=parametername; + end end + rethrow(me2); end Index: /issm/trunk/src/m/classes/public/plot/applyoptions.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/applyoptions.m (revision 26) +++ /issm/trunk/src/m/classes/public/plot/applyoptions.m (revision 27) @@ -44,4 +44,10 @@ if ~isnan(options_structure.view), view(options_structure.view); +else + if strcmpi(md.type,'3d') & isnan(options_structure.layer), + view(3); + else + view(2); + end end @@ -162,4 +168,9 @@ end +%streamliness +if iscell(options_structure.streamlines) | ~isnan(options_structure.streamlines), + plot_streamlines(md,options_structure); +end + %contours if iscell(options_structure.contourlevels) | ~isnan(options_structure.contourlevels), Index: /issm/trunk/src/m/classes/public/plot/imagescnan.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/imagescnan.m (revision 26) +++ /issm/trunk/src/m/classes/public/plot/imagescnan.m (revision 27) @@ -80,5 +80,5 @@ % Copyright 2008 Carlos Adrian Vargas Aguilera -% $Revision: 1.1 $ $Date: 2009/03/24 21:05:20 $ +% $Revision: 1.1 $ $Date: 2009/04/03 22:56:05 $ % Written by Index: /issm/trunk/src/m/classes/public/plot/parse_options.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/parse_options.m (revision 26) +++ /issm/trunk/src/m/classes/public/plot/parse_options.m (revision 27) @@ -17,4 +17,12 @@ end +%density +densityvalues=findarg(optionstring,'density'); +if ~isempty(densityvalues), + options_struct.density=abs(ceil(densityvalues(1).value)); +else + options_struct.density=NaN; +end + %Section profile sectionvalues=findarg(optionstring,'sectionvalue'); @@ -38,10 +46,10 @@ %Show section -showsections=findarg(optionstring,'showsection'); -if ~isempty(showsections), - if ischar(showsections(1).value), - options_struct.showsection=showsections(1).value; - else - options_struct.showsection=0; +showsection=findarg(optionstring,'showsection'); +if ~isempty(showsection), + if ischar(showsection(1).value) & strcmpi(showsection(1).value,'on'), + options_struct.showsection=4; + else + options_struct.showsection=showsection(1).value; end else @@ -115,4 +123,12 @@ else options_struct.smooth=NaN; +end + +%streamlines +streamlines_values=findarg(optionstring,'streamlines'); +if ~isempty(streamlines_values), + options_struct.streamlines=streamlines_values.value; +else + options_struct.streamlines=NaN; end @@ -226,9 +242,5 @@ options_struct.view=viewvalues.value; else - if strcmpi(md.type,'3d') & isnan(options_struct.layer), - options_struct.view=3; - else - options_struct.view=2; - end + options_struct.view=NaN; end Index: /issm/trunk/src/m/classes/public/plot/plot_basaldrag.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plot_basaldrag.m (revision 26) +++ /issm/trunk/src/m/classes/public/plot/plot_basaldrag.m (revision 27) @@ -1,7 +1,7 @@ -function plot_basaldrag(md,options_structure,width,i); +function plot_basaldrag(md,options_structure,width,i,type); %PLOT_BASALDRAG - plot basal drag % % Usage: -% plot_basaldrag(md,options_structure,width,i); +% plot_basaldrag(md,options_structure,width,i,type); % % See also: PLOTMODEL @@ -20,32 +20,37 @@ %compute horizontal velocity -ub=sqrt(md.vx.^2+md.vy.^2)/md.yts; +if strcmpi(type,'basal_drag') + ub=sqrt(md.vx.^2+md.vy.^2)/md.yts; +elseif strcmpi(type,'basal_dragx') + ub=md.vx/md.yts; +elseif strcmpi(type,'basal_dragy') + ub=md.vy/md.yts; +end %compute basal drag drag=(md.g*(md.rho_ice*md.thickness+md.rho_water*md.bed)).^r.*(md.drag).^2.*ub.^s/1000; -%process data and model -[x y z elements is2d]=processmesh(md,options_structure); -[basal_drag isongrid]=processdata(md,drag,options_structure); +%Figure out if this is a Section plot +if ~isnan(options_structure.sectionvalue) + plot_section(md,drag,options_structure,width,i); + return; +else -%plot mesh quivervel -subplot(width,width,i); + %process data and model + [x y z elements is2d]=processmesh(md,options_structure); + [basal_drag isongrid]=processdata(md,drag,options_structure); -%edgecolor? -if ~isnan(options_structure.edgecolor), - edgecolor=options_structure.edgecolor; -else - edgecolor='none'; + %plot basaldrag + subplot(width,width,i); + plot_unit(x,y,z,elements,basal_drag,isongrid,is2d,options_structure); + + %apply options + if isnan(options_structure.title) + options_structure.title='Basal drag [kPa]'; + end + if isnan(options_structure.view) + options_structure.view=2; + end + applyoptions(md,basal_drag,options_structure); + end - -%plot basal basal drag -A=elements(:,1); B=elements(:,2); C=elements(:,3); - -patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData',basal_drag,'FaceColor','interp','EdgeColor',edgecolor); - -%apply options -if isnan(options_structure.title) - options_structure.title='Basal drag [kPa]'; -end -options_structure.view=2; -applyoptions(md,basal_drag,options_structure); Index: /issm/trunk/src/m/classes/public/plot/plot_boundaries.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plot_boundaries.m (revision 26) +++ /issm/trunk/src/m/classes/public/plot/plot_boundaries.m (revision 27) @@ -41,3 +41,6 @@ options_structure.colorbar=0; end +if isnan(options_structure.view) + options_structure.view=2; +end applyoptions(md,[],options_structure); Index: /issm/trunk/src/m/classes/public/plot/plot_contour.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plot_contour.m (revision 26) +++ /issm/trunk/src/m/classes/public/plot/plot_contour.m (revision 27) @@ -46,5 +46,5 @@ %initialization of some variables -numberofelements=md.numberofelements; +numberofelements=size(index,1); elementslist=1:numberofelements; c=[]; @@ -231,5 +231,5 @@ %labels? -if ~strcmpi(options_structure.contourticks,'off') +if (~strcmpi(options_structure.contourticks,'off') & ~isempty(c) & ~isempty(h)) if ~isnan(options_structure.contourcolor) clabel(c,h,'color',color); Index: /issm/trunk/src/m/classes/public/plot/plot_dealer.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plot_dealer.m (revision 27) +++ /issm/trunk/src/m/classes/public/plot/plot_dealer.m (revision 27) @@ -0,0 +1,129 @@ +function plot_dealer(md,optionstring,width,i); +%PLOT_DEALER - distribute the plots, called by plotmodel +% +% Usage: +% plot_dealer(md,optionstring,width,i); +% +% See also: PLOTMODEL, PLOT_UNIT + +%parse options and get a structure of options. +options_structure=parse_options(md,optionstring); + +%get data to be displayed +data=findarg(optionstring,'data');data=data.value; + +%figure out if this is a special plot +if ischar(data), + + switch data, + + case 'boundaries', + plot_boundaries(md,options_structure,width,i); + return; + case 'elementnumbering', + plot_elementnumbering(md,options_structure,width,i); + return; + case 'highlightelements', + plot_highlightelements(md,options_structure,width,i); + return; + case 'segmentnumbering', + plot_segmentnumbering(md,options_structure,width,i); + return; + case 'importancefactors', + plot_importancefactors(md,options_structure,width,i); + return; + case 'elements_type', + plot_elementstype(md,options_structure,width,i); + return; + case 'gridnumbering', + plot_gridnumbering(md,options_structure,width,i); + return; + case 'highlightgrids', + plot_highlightgrids(md,options_structure,width,i); + return; + case {'basal_drag','basal_dragx','basal_dragy'}, + plot_basaldrag(md,options_structure,width,i,data); + return; + case 'driving_stress', + plot_drivingstress(md,options_structure,width,i); + return; + case 'mesh', + plot_mesh(md,options_structure,width,i); + return; + case 'penalties', + plot_penalties(md,options_structure,width,i); + return; + case 'quiver', + plot_quiver(md,options_structure,width,i); + return; + case 'quiver3', + plot_quiver3(md,options_structure,width,i); + return; + case 'quivervel', + plot_quivervel(md,options_structure,width,i); + return; + case 'riftvel', + plot_riftvel(md,options_structure,width,i); + return; + case 'riftrelvel', + plot_riftrelvel(md,options_structure,width,i); + return; + case 'riftpenetration', + plot_riftpenetration(md,options_structure,width,i); + return; + case 'sarpwr', + plot_sarpwr(md,options_structure,width,i) + return + case {'segmentonneumann_diag','segmentonneumann_prog'} + plot_segmentonneumann(md,options_structure,width,i,data) + return + case {'strainrate_tensor','strainrate','strainrate_principal','strainrate_principalaxis1','strainrate_principalaxis2','strainrate_principalaxis3',... + 'stress_tensor','stress','stress_principal','stress_principalaxis1','stress_principalaxis2','stress_principalaxis3',... + 'deviatoricstress_tensor','deviatoricstress','deviatoricstress_principal','deviatoricstress_principalaxis1','deviatoricstress_principalaxis2','deviatoricstress_principalaxis3'}, + plot_tensor(md,options_structure,width,i,data); + return; + case 'thermaltransient_results', + plot_thermaltransient_results(md,options_structure,width,i); + return; + case 'transient_movie', + plot_transient_movie(md,options_structure,width,i); + return; + case 'transient_results', + plot_transient_results(md,options_structure,width,i); + return; + + otherwise, + + if isfield(struct(md),data) + data=eval(['md.' data ';']); + else + error('plot error message: data provided not supported yet. Type plotdoc for help'); + end + end +end + +%Figure out if this is a semi-transparent plot. +if ~isnan(options_structure.overlay), + plot_overlay(md,data,options_structure,width,i); + return; +end + +%Figure out if this is a Section plot +if ~isnan(options_structure.sectionvalue) + plot_section(md,data,options_structure,width,i); + return; +end + +%process data and model +[x y z elements is2d]=processmesh(md,options_structure); +[data isongrid]=processdata(md,data,options_structure); + +%standard plot: +subplot(width,width,i); +plot_unit(x,y,z,elements,data,isongrid,is2d,options_structure); + +%apply all options +if isnan(options_structure.shading) & isnan(options_structure.edgecolor) & size(data,1)==md.numberofgrids, + options_structure.shading='interp'; +end +applyoptions(md,data,options_structure); Index: /issm/trunk/src/m/classes/public/plot/plot_drivingstress.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plot_drivingstress.m (revision 26) +++ /issm/trunk/src/m/classes/public/plot/plot_drivingstress.m (revision 27) @@ -5,5 +5,5 @@ % plot_drivingstress(md,options_structure,width,i); % -% See also: PLOTMODEL +% See also: PLOTMODEL, PLOT_UNIT, PLOT_DEALER %get driving stress @@ -17,44 +17,5 @@ %plot mesh quivervel subplot(width,width,i); - -%edgecolor? -if ~isnan(options_structure.edgecolor), - edgecolor=options_structure.edgecolor; -else - edgecolor='none'; -end - -%element data -if ~isongrid - if is2d - A=elements(:,1); B=elements(:,2); C=elements(:,3); - patch( 'Faces', [A B C], 'Vertices', [x y z],'CData', dstress,'FaceColor','flat','EdgeColor',edgecolor); - else - A=elements(:,1); B=elements(:,2); C=elements(:,3); D=elements(:,4); E=elements(:,5); F=elements(:,6); - patch( 'Faces', [A B C], 'Vertices', [x y z],'CData', dstress,'FaceColor','flat','EdgeColor',edgecolor); - patch( 'Faces', [D E F], 'Vertices', [x y z],'CData', dstress,'FaceColor','flat','EdgeColor',edgecolor); - patch( 'Faces', [A B E D], 'Vertices', [x y z],'CData', dstress,'FaceColor','flat','EdgeColor',edgecolor); - patch( 'Faces', [B E F C ], 'Vertices', [x y z],'CData', dstress,'FaceColor','flat','EdgeColor',edgecolor); - patch( 'Faces', [C A D F ], 'Vertices', [x y z],'CData', dstress,'FaceColor','flat','EdgeColor',edgecolor); - end -%grid data -elseif isongrid - if is2d - A=elements(:,1); B=elements(:,2); C=elements(:,3); - patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', dstress,'FaceColor','interp','EdgeColor',edgecolor); - else - if options_structure.layer>=1, - A=elements(:,1); B=elements(:,2); C=elements(:,3); - patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', dstress,'FaceColor','interp','EdgeColor',edgecolor); - else - A=elements(:,1); B=elements(:,2); C=elements(:,3); D=elements(:,4); E=elements(:,5); F=elements(:,6); - patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', dstress,'FaceColor','interp','EdgeColor',edgecolor); - patch( 'Faces', [D E F], 'Vertices', [x y z],'FaceVertexCData', dstress,'FaceColor','interp','EdgeColor',edgecolor); - patch( 'Faces', [A B E D], 'Vertices', [x y z],'FaceVertexCData', dstress,'FaceColor','interp','EdgeColor',edgecolor); - patch( 'Faces', [B E F C ], 'Vertices', [x y z],'FaceVertexCData', dstress,'FaceColor','interp','EdgeColor',edgecolor); - patch( 'Faces', [C A D F ], 'Vertices', [x y z],'FaceVertexCData', dstress,'FaceColor','interp','EdgeColor',edgecolor); - end - end -end +plot_unit(x,y,z,elements,dstress,isongrid,is2d,options_structure) %apply options @@ -62,3 +23,6 @@ options_structure.title='Driving stress [kPa]'; end +if isnan(options_structure.view) + options_structure.view=2; +end applyoptions(md,dstress,options_structure); Index: /issm/trunk/src/m/classes/public/plot/plot_elementnumbering.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plot_elementnumbering.m (revision 26) +++ /issm/trunk/src/m/classes/public/plot/plot_elementnumbering.m (revision 27) @@ -5,5 +5,5 @@ % plot_elementnumbering(md,options_structure,width,i); % -% See also: PLOTMODEL +% See also: PLOTMODEL, PLOT_UNIT, PLOT_DEALER subplot(width,width,i); Index: /issm/trunk/src/m/classes/public/plot/plot_elementstype.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plot_elementstype.m (revision 26) +++ /issm/trunk/src/m/classes/public/plot/plot_elementstype.m (revision 27) @@ -76,10 +76,10 @@ if ~isempty(posS) A=elements(posS,1); B=elements(posS,2); C=elements(posS,3); D=elements(posS,4); E=elements(posS,5); F=elements(posS,6); - % p4=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData', stokesenum,'FaceColor','flat','EdgeColor',edgecolor,'EdgeAlpha',alpha,'FaceAlpha',alpha); - % patch( 'Faces', [D E F], 'Vertices', [x y z],'CData', stokesenum,'FaceColor','flat','EdgeColor',edgecolor,'EdgeAlpha',alpha,'FaceAlpha',alpha); - % patch( 'Faces', [A B E D], 'Vertices', [x y z],'CData', stokesenum,'FaceColor','flat','EdgeColor',edgecolor,'EdgeAlpha',alpha,'FaceAlpha',alpha); - % patch( 'Faces', [B E F C ], 'Vertices', [x y z],'CData', stokesenum,'FaceColor','flat','EdgeColor',edgecolor,'EdgeAlpha',alpha,'FaceAlpha',alpha); - % patch( 'Faces', [C A D F ], 'Vertices', [x y z],'CData', stokesenum,'FaceColor','flat','EdgeColor',edgecolor,'EdgeAlpha',alpha,'FaceAlpha',alpha); - % legend([p1 p2 p3 p4],'Hutter''s elements','MacAyeal''s elements','Pattyn''s elements','Stokes''s elements'); + p4=patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceColor','none','EdgeColor','black'); + patch( 'Faces', [D E F], 'Vertices', [x y z],'FaceColor','none','EdgeColor','black'); + patch( 'Faces', [A B E D], 'Vertices', [x y z],'FaceColor','none','EdgeColor','black'); + patch( 'Faces', [B E F C ], 'Vertices', [x y z],'FaceColor','none','EdgeColor','black'); + patch( 'Faces', [C A D F ], 'Vertices', [x y z],'FaceColor','none','EdgeColor','black'); + legend([p1 p2 p3 p4],'Hutter''s elements','MacAyeal''s elements','Pattyn''s elements','Stokes''s elements'); else legend([p1 p2 p3],'Hutter''s elements','MacAyeal''s elements','Pattyn''s elements'); Index: /issm/trunk/src/m/classes/public/plot/plot_penalties.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plot_penalties.m (revision 26) +++ /issm/trunk/src/m/classes/public/plot/plot_penalties.m (revision 27) @@ -22,5 +22,5 @@ if ~strcmpi(md.type,'3d'), error('no penalties to plot for ''2d'' model'); -elseif isnan(md.penalties) | isempty(md.penalties) +elseif isempty(md.penalties), disp('no penalty applied in this model'); return; Index: /issm/trunk/src/m/classes/public/plot/plot_quiver.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plot_quiver.m (revision 26) +++ /issm/trunk/src/m/classes/public/plot/plot_quiver.m (revision 27) @@ -16,4 +16,10 @@ %quiver 2d +if ~isnan(options_structure.density) + x=x(1:options_structure.density:end); + y=y(1:options_structure.density:end); + vx=vx(1:options_structure.density:end); + vy=vy(1:options_structure.density:end); +end quiver(x,y,vx,vy); Index: /issm/trunk/src/m/classes/public/plot/plot_quiver3.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plot_quiver3.m (revision 26) +++ /issm/trunk/src/m/classes/public/plot/plot_quiver3.m (revision 27) @@ -17,4 +17,12 @@ %quiver 3d +if ~isnan(options_structure.density) + x=x(1:options_structure.density:end); + y=y(1:options_structure.density:end); + z=z(1:options_structure.density:end); + vx=vx(1:options_structure.density:end); + vy=vy(1:options_structure.density:end); + vz=vz(1:options_structure.density:end); +end quiver3(x,y,z,vx,vy,vz); Index: /issm/trunk/src/m/classes/public/plot/plot_quivervel.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plot_quivervel.m (revision 26) +++ /issm/trunk/src/m/classes/public/plot/plot_quivervel.m (revision 27) @@ -5,5 +5,5 @@ % plot_quivervel(md,options_structure,width,i); % -% See also: PLOTMODEL +% See also: PLOTMODEL, PLOT_UNIT, PLOT_DEALER %process data and model @@ -15,19 +15,21 @@ end -%edgecolor? -if ~isnan(options_structure.edgecolor), - edgecolor=options_structure.edgecolor; -else - edgecolor='none'; +%density +if ~isnan(options_structure.density) + x=x(1:options_structure.density:end); + y=y(1:options_structure.density:end); + z=z(1:options_structure.density:end); + vx=vx(1:options_structure.density:end); + vy=vy(1:options_structure.density:end); + vz=vz(1:options_structure.density:end); end %plot subplot(width,width,i); +colormap('default') +plot_unit(x,y,z,elements,sqrt(vx.^2+vy.^2),isongrid,is2d,options_structure) if is2d hold on - colormap('default') - A=elements(:,1); B=elements(:,2); C=elements(:,3); - patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor); colorbar; quiver(x,y,vx,vy,'k') @@ -35,11 +37,4 @@ else hold on - colormap('default') - A=elements(:,1); B=elements(:,2); C=elements(:,3); D=elements(:,4); E=elements(:,5); F=elements(:,6); - patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor); - patch( 'Faces', [D E F], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor); - patch( 'Faces', [A B E D], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor); - patch( 'Faces', [B E F C ], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor); - patch( 'Faces', [C A D F ], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor); colorbar; quiver3(x,y,z,vx,vy,vz); Index: /issm/trunk/src/m/classes/public/plot/plot_section.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plot_section.m (revision 26) +++ /issm/trunk/src/m/classes/public/plot/plot_section.m (revision 27) @@ -8,5 +8,5 @@ %How many subplots? -if strcmpi(options_structure.showsection,'yes') +if ~isnan(options_structure.showsection) %Compute the indexes of the 2 plots (one for the sectionvalue and one for showsection @@ -20,34 +20,12 @@ end -%check on arguments -if (iscell(data) | isempty(data)), - error('plot error message: data provided is empty'); -end -if md.numberofgrids==size(md.elements,1), - error('plot error message: the number of elements is the same as the number of grids! cannot plot anything with model/plot, use matlab/plot instead') -end -%smoothing? -if strcmpi(options_structure.smooth,'yes') & length(data)==md.numberofelements - data=elementstogrids(md,data); -end +%process data and model +[x_m y_m z_m elements_m is2d]=processmesh(md,options_structure); +[data isongrid]=processdata(md,data,options_structure); -%layer projection? -if ~isnan(options_structure.layer) & options_structure.layer>=1, - data=project2d(md,data,options_structure.layer); %project onto 2d mesh - %we modify the mesh temporarily to a 2d mesh from which the 3d mesh was extruded. - md.x=md.x2d; - md.y=md.y2d; - md.z=md.z2d; - md.elements=md.elements2d; - md.elements_type=md.elements_type2d; - md.type='2d'; -end - -%units -if ~isnan(options_structure.unitmultiplier), - md.x=md.x*options_structure.unitmultiplier; - md.y=md.y*options_structure.unitmultiplier; - md.z=md.z*options_structure.unitmultiplier; +%replug x and y onto model so that SectionValue treats the problem correctly +if ~isnan(options_structure.layer) + md.x=md.x2d; md.y=md.y2d; md.elements=md.elements2d; md.type='2d'; end @@ -61,10 +39,39 @@ %Compute section value -[elements,x,y,z,s,data]=SectionValues(md,data,options_structure.sectionvalue,resolution); +[elements,x,y,z,s,data_s]=SectionValues(md,data,options_structure.sectionvalue,resolution); -if strcmpi(md.type,'2d') +if is2d %plot section value subplot(width,width,index1) - plot(s,data) + plot(s,data_s) + + %Show Section if requested by user + if ~isnan(options_structure.showsection) + + %compute number of labels + numlabels=min(options_structure.showsection,length(s)); + shift=fix(length(s)/numlabels); + + %plot labels on current graph + hold on + text(s(1),data_s(1),'1','backgroundcolor',[0.8 0.9 0.8]) + for i=2:numlabels-1 + text(s(1+(i-1)*shift),data_s(1+(i-1)*shift),num2str(i),'backgroundcolor',[0.8 0.9 0.8]) + end + text(s(end),data_s(end),'end','backgroundcolor',[0.8 0.9 0.8]) + + %plot section only with labels + subplot(width,width,index2) + plot_unit(x_m,y_m,z_m,elements_m,data,isongrid,is2d,options_structure) + hold on + text(x(1),y(1),'1','backgroundcolor',[0.8 0.9 0.8]) + for i=2:numlabels-1 + text(x(1+(i-1)*shift),y(1+(i-1)*shift),num2str(i),'backgroundcolor',[0.8 0.9 0.8]) + end + text(x(end),y(end),'end','backgroundcolor',[0.8 0.9 0.8]) + plot(x,y,'-r') + axis([min(md.x)-1 max(md.x)+1 min(md.y)-1 max(md.y)+1]) + view(2) + end else %plot section value @@ -73,10 +80,68 @@ subplot(width,width,index1) A=elements(:,1); B=elements(:,2); C=elements(:,3); D=elements(:,4); - patch( 'Faces', [A B C D], 'Vertices', [s z zeros(length(s),1)],'FaceVertexCData',data,'FaceColor','interp','EdgeColor','none'); + patch( 'Faces', [A B C D], 'Vertices', [s z zeros(length(s),1)],'FaceVertexCData',data_s,'FaceColor','interp','EdgeColor','none'); + + %Show Section if requested by user + if ~isnan(options_structure.showsection) + + %compute number of labels + numlabels=min(options_structure.showsection,length(s)); + shift=fix(length(s)/numlabels); + + %plot labels on current graph + hold on + text(s(1),z(1),'1','backgroundcolor',[0.8 0.9 0.8]) + for i=2:numlabels-1 + text(s(1+(i-1)*shift),z(1+(i-1)*shift),num2str(i),'backgroundcolor',[0.8 0.9 0.8]) + end + text(s(end),z(end),'end','backgroundcolor',[0.8 0.9 0.8]) + + %plot section only with labels + subplot(width,width,index2) + plot_unit(x_m,y_m,z_m,elements_m,data,isongrid,is2d,options_structure) + hold on + text(x(1),y(1),'1','backgroundcolor',[0.8 0.9 0.8]) + for i=2:numlabels-1 + text(x(1+(i-1)*shift),y(1+(i-1)*shift),num2str(i),'backgroundcolor',[0.8 0.9 0.8]) + end + text(x(end),y(end),'end','backgroundcolor',[0.8 0.9 0.8]) + plot(x,y,'-r') + axis([min(md.x)-1 max(md.x)+1 min(md.y)-1 max(md.y)+1]) + view(2) + end else subplot(width,width,index1) A=elements(:,1); B=elements(:,2); C=elements(:,3); D=elements(:,4); - patch( 'Faces', [A B C D], 'Vertices', [x y z],'FaceVertexCData',data,'FaceColor','interp','EdgeColor','none'); + patch( 'Faces', [A B C D], 'Vertices', [x y z],'FaceVertexCData',data_s,'FaceColor','interp','EdgeColor','none'); view(3) + + %Show Section if requested by user + if ~isnan(options_structure.showsection) + + %compute number of labels + numlabels=min(options_structure.showsection,length(s)); + shift=fix(length(x)/numlabels); + + %plot labels on current graph + hold on + text(x(1),y(1),z(1),'1','backgroundcolor',[0.8 0.9 0.8]) + for i=2:numlabels-1 + text(x(1+(i-1)*shift),y(1+(i-1)*shift),z(1+(i-1)*shift),num2str(i),'backgroundcolor',[0.8 0.9 0.8]) + end + text(x(end),y(end),z(end),'end','backgroundcolor',[0.8 0.9 0.8]) + + %plot section only with labels + subplot(width,width,index2) + plot_unit(x_m,y_m,z_m,elements_m,data,isongrid,is2d,options_structure) + hold on + text(x(1),y(1),'1','backgroundcolor',[0.8 0.9 0.8]) + for i=2:numlabels-1 + text(x(1+(i-1)*shift),y(1+(i-1)*shift),num2str(i),'backgroundcolor',[0.8 0.9 0.8]) + end + text(x(end),y(end),'end','backgroundcolor',[0.8 0.9 0.8]) + plot(x,y,'-r') + axis([min(md.x)-1 max(md.x)+1 min(md.y)-1 max(md.y)+1]) + view(2) + end end end @@ -96,13 +161,2 @@ end applyoptions(md,[],options_structure); - -%plot section if requested by user -if strcmpi(options_structure.showsection,'yes') - subplot(width,width,index2) - hold on - text(x(1),y(1),'1','backgroundcolor',[0.8 0.9 0.8]) - text(x(end),y(end),'end','backgroundcolor',[0.8 0.9 0.8]) - plot(x,y) - axis([min(md.x) max(md.x) min(md.y) max(md.y)]) - view(2) -end Index: /issm/trunk/src/m/classes/public/plot/plot_streamlines.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plot_streamlines.m (revision 27) +++ /issm/trunk/src/m/classes/public/plot/plot_streamlines.m (revision 27) @@ -0,0 +1,156 @@ +function plot_streamlines(md,options_structure) +%PLOT_STREAMLINES - plot stream lines on a figure +% +% Usage: +% plot_streamlines(md,options_structure) + +%some options +maxiter=200; %maximum number of iterations +precision=1; %division of each segment (higer precision increases number of segments) + +%process data and model +[x y z index is2d]=processmesh(md,options_structure); +[u isongrid]=processdata(md,md.vx,options_structure); +[v isongrid]=processdata(md,md.vy,options_structure); + +%some checks +if ~is2d, + disp('plot_streamlines error: streamlines option not supported for 3d plots. Project on a layer') + return +end + +%take velocity for each element +u=u(index)*[1;1;1]/3; +v=v(index)*[1;1;1]/3; + +%initialize flowpath +if iscell(options_structure.streamlines) + x0=[]; y0=[]; + for i=1:size(options_structure.streamlines,2) + coord=options_structure.streamlines{i}; + x0=[x0;coord(1)]; y0=[y0;coord(2)]; + end +else + x0=x(1:ceil(length(x)/options_structure.streamlines):end); + y0=y(1:ceil(length(x)/options_structure.streamlines):end); +end + +%check seed points +tria=tsearch(x,y,index,x0,y0); +pos=find(isnan(tria)); +%seems to be a bug here: must do the test several times... +while ~isempty(pos) + x0(pos)=[]; + y0(pos)=[]; + tria=tsearch(x,y,index,x0,y0); + pos=find(isnan(tria)); +end + +%initialize other variables +N=length(x0); +X=x0; Y=y0; +flowpath=struct('xc',cell(N,1),'yc',cell(N,1),'done',cell(N,1)); +for i=1:N, + flowpath(i).xc=x0(i); + flowpath(i).yc=y0(i); +end +done=zeros(N,1); + +%get avegared length of each element +length_tria=1/3*(sqrt( (x(index(:,1))-x(index(:,2))).^2+(y(index(:,1))-y(index(:,2))).^2 )+... + sqrt((x(index(:,1))-x(index(:,3))).^2+(y(index(:,1))-y(index(:,3))).^2 )+... + sqrt((x(index(:,2))-x(index(:,3))).^2+(y(index(:,2))-y(index(:,3))).^2 )); + +%initialization: +counter=1; + +while any(~done) + + %find current triangle + queue=find(~done); + tria=tsearch(x,y,index,X(queue),Y(queue)); + + %check that the point is actually inside a triangle of the mesh + listnan=find(isnan(tria)); + for i=1:length(listnan) + %remove the last point + flowpath(queue(listnan(i))).xc(end)=[]; + flowpath(queue(listnan(i))).yc(end)=[]; + done(queue(listnan(i)))=1; + end + tria(listnan)=[]; + queue(listnan)=[]; + + %velocity of the current triangle and norm it + ut=u(tria); vt=v(tria); normv=sqrt(ut.^2+vt.^2); + ut=ut./normv;vt=vt./normv; + + %check counter + if counter>maxiter + disp(['Maximum number of iterations (' num2str(maxiter) ') reached while going forward']) + break + end + counter=counter+1; + + %remove stagnant point + done(queue(find(ut==0 & vt==0)))=1; + + %build next point + for i=1:length(queue) + X(queue(i))=flowpath(queue(i)).xc(end)+ut(i)*length_tria(tria(i))/precision; + Y(queue(i))=flowpath(queue(i)).yc(end)+vt(i)*length_tria(tria(i))/precision; + flowpath(queue(i)).xc=[flowpath(queue(i)).xc;flowpath(queue(i)).xc(end)+ut(i)*length_tria(tria(i))/precision]; + flowpath(queue(i)).yc=[flowpath(queue(i)).yc;flowpath(queue(i)).yc(end)+vt(i)*length_tria(tria(i))/precision]; + end +end + +%same process but reverse (vel=-vel) to have a vcomplete flow line +counter=1; +X=x0; Y=y0; +done=zeros(N,1); + +while any(~done) + + %find current triangle + queue=find(~done); + tria=tsearch(x,y,index,X(queue),Y(queue)); + + %check that the point is actually inside a triangle of the mesh + listnan=find(isnan(tria)); + for i=1:length(listnan) + %remove the last point + flowpath(queue(listnan(i))).xc(1)=[]; + flowpath(queue(listnan(i))).yc(1)=[]; + done(queue(listnan(i)))=1; + end + tria(listnan)=[]; + queue(listnan)=[]; + + %velocity of the current triangle and norm it + ut=-u(tria); vt=-v(tria); normv=sqrt(ut.^2+vt.^2); + ut=ut./normv;vt=vt./normv; + + %check counter + if counter>maxiter + disp(['Maximum number of iterations (' num2str(maxiter) ') reached while going backward']) + break + end + counter=counter+1; + + %remove stagnant point + done(queue(find(ut==0 & vt==0)))=1; + + %build next point + for i=1:length(queue) + X(queue(i))=flowpath(queue(i)).xc(1)+ut(i)*length_tria(tria(i))/precision; + Y(queue(i))=flowpath(queue(i)).yc(1)+vt(i)*length_tria(tria(i))/precision; + flowpath(queue(i)).xc=[flowpath(queue(i)).xc(1)+ut(i)*length_tria(tria(i))/precision; flowpath(queue(i)).xc]; + flowpath(queue(i)).yc=[flowpath(queue(i)).yc(1)+vt(i)*length_tria(tria(i))/precision; flowpath(queue(i)).yc]; + end +end + +%plot +hold on +for i=1:length(flowpath) + patch('Xdata',[flowpath(i).xc;NaN],'Ydata',[flowpath(i).yc;NaN],'facecolor','none','edgecolor','y'); +end Index: /issm/trunk/src/m/classes/public/plot/plot_tensor_components.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plot_tensor_components.m (revision 26) +++ /issm/trunk/src/m/classes/public/plot/plot_tensor_components.m (revision 27) @@ -5,298 +5,67 @@ % plot_tensor_components(md,options_structure,width,i,tensor,type,plot_option); % -% See also: PLOTMODEL +% See also: PLOTMODEL, PLOT_UNIT, PLOT_DEALER - %Compute the indexes of the components plots - upperplots=fix((i-1)/width); - if upperplots==0, leftplots=i-1; else leftplots=i-width*upperplots-1; end - if strcmpi(md.type,'2d')%3 components -> 3 indexes - index1=4*width*upperplots+2*leftplots+1; - index2=index1+1; - index3=index1+width*2; - elseif strcmpi(md.type,'3d')%6 components -> 6 indexes - index1=3*3*width*upperplots+3*leftplots+1; - index2=index1+1; - index3=index1+2; - index4=index1+width*3; - index5=index4+1; - index6=index4+2; - end - - %smoothing? - if strcmpi(options_structure.smooth,'yes') & length(tensor.xx)==md.numberofelements - tensor.xx=elementstogrids(md,tensor.xx); - tensor.yy=elementstogrids(md,tensor.yy); - tensor.xy=elementstogrids(md,tensor.xy); - if strcmpi(md.type,'3d') - tensor.zz=elementstogrids(md,tensor.zz); - tensor.yz=elementstogrids(md,tensor.yz); - tensor.xz=elementstogrids(md,tensor.xz); - end - end - - - %layer projection? - if ~isnan(options_structure.layer) & options_structure.layer>=1, - tensor.xx=project2d(md,tensor.xx,options_structure.layer); %project onto 2d mesh - tensor.yy=project2d(md,tensor.yy,options_structure.layer); %project onto 2d mesh - tensor.zz=project2d(md,tensor.zz,options_structure.layer); %project onto 2d mesh - tensor.xy=project2d(md,tensor.xy,options_structure.layer); %project onto 2d mesh - tensor.xz=project2d(md,tensor.xz,options_structure.layer); %project onto 2d mesh - tensor.yz=project2d(md,tensor.yz,options_structure.layer); %project onto 2d mesh - %we modify the mesh temporarily to a 2d mesh from which the 3d mesh was extruded. - md.x=md.x2d; - md.y=md.y2d; - md.z=md.z2d; - md.elements=md.elements2d; - md.elements_type=md.elements_type2d; - end - - %units - if ~isnan(options_structure.unitmultiplier), - md.x=md.x*options_structure.unitmultiplier; - md.y=md.y*options_structure.unitmultiplier; - md.z=md.z*options_structure.unitmultiplier; - end - - if length(tensor.xx)==length(md.elements), - - if ~isnan(options_structure.noiceshelf) & options_structure.noiceshelf, - tensor.xx(find(md.elementoniceshelf))=NaN; - tensor.yy(find(md.elementoniceshelf))=NaN; - tensor.xy(find(md.elementoniceshelf))=NaN; - if strcmpi(md.type,'3d') - tensor.zz(find(md.elementoniceshelf))=NaN; - tensor.yz(find(md.elementoniceshelf))=NaN; - tensor.xz(find(md.elementoniceshelf))=NaN; - end - end - if ~isnan(options_structure.noicesheet) & options_structure.noicesheet, - tensor.xx(find(~md.elementoniceshelf))=NaN; - tensor.yy(find(~md.elementoniceshelf))=NaN; - tensor.xy(find(~md.elementoniceshelf))=NaN; - if strcmpi(md.type,'3d') - tensor.zz(find(~md.elementoniceshelf))=NaN; - tensor.yz(find(~md.elementoniceshelf))=NaN; - tensor.xz(find(~md.elementoniceshelf))=NaN; - end - end - - if (strcmpi(md.type,'2d')), - A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3); - - subplot(2*width,2*width,index1) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.xx,'FaceColor','flat','EdgeColor','black'); - Apply_options_tensor(options_structure,type,'xx') - - subplot(2*width,2*width,index2) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.yy,'FaceColor','flat','EdgeColor','black'); - Apply_options_tensor(options_structure,type,'yy') - - subplot(2*width,2*width,index3) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.xy,'FaceColor','flat','EdgeColor','black'); - Apply_options_tensor(options_structure,type,'xy') - else - if options_structure.layer>=1, - A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3); - - subplot(3*width,3*width,index1) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData',tensor.xx ,'FaceColor','flat','EdgeColor','black'); - Apply_options_tensor(options_structure,type,'xx') - - subplot(3*width,3*width,index2) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData',tensor.yy ,'FaceColor','flat','EdgeColor','black'); - Apply_options_tensor(options_structure,type,'yy') - - subplot(3*width,3*width,index3) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData',tensor.zz ,'FaceColor','flat','EdgeColor','black'); - Apply_options_tensor(options_structure,type,'zz') - - subplot(3*width,3*width,index4) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData',tensor.xy ,'FaceColor','flat','EdgeColor','black'); - Apply_options_tensor(options_structure,type,'xy') - - subplot(3*width,3*width,index5) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData',tensor.xz ,'FaceColor','flat','EdgeColor','black'); - Apply_options_tensor(options_structure,type,'xz') - - subplot(3*width,3*width,index6) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData',tensor.yz ,'FaceColor','flat','EdgeColor','black'); - Apply_options_tensor(options_structure,type,'yz') - - else - A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3); D=md.elements(:,4); E=md.elements(:,5); F=md.elements(:,6); - - subplot(3*width,3*width,index1) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.xx,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'CData', tensor.xx,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'CData', tensor.xx,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'CData', tensor.xx,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'CData', tensor.xx,'FaceColor','flat','EdgeColor','black'); - Apply_options_tensor(options_structure,type,'xx') - - subplot(3*width,3*width,index2) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.yy,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'CData', tensor.yy,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'CData', tensor.yy,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'CData', tensor.yy,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'CData', tensor.yy,'FaceColor','flat','EdgeColor','black'); - Apply_options_tensor(options_structure,type,'yy') - - subplot(3*width,3*width,index3) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.zz,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'CData', tensor.zz,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'CData', tensor.zz,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'CData', tensor.zz,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'CData', tensor.zz,'FaceColor','flat','EdgeColor','black'); - Apply_options_tensor(options_structure,type,'zz') - - subplot(3*width,3*width,index4) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.xy,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'CData', tensor.xy,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'CData', tensor.xy,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'CData', tensor.xy,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'CData', tensor.xy,'FaceColor','flat','EdgeColor','black'); - Apply_options_tensor(options_structure,type,'xy') - - subplot(3*width,3*width,index5) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.xz,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'CData', tensor.xz,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'CData', tensor.xz,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'CData', tensor.xz,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'CData', tensor.xz,'FaceColor','flat','EdgeColor','black'); - Apply_options_tensor(options_structure,type,'xz') - - subplot(3*width,3*width,index6) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.yz,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'CData', tensor.yz,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'CData', tensor.yz,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'CData', tensor.yz,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'CData', tensor.yz,'FaceColor','flat','EdgeColor','black'); - Apply_options_tensor(options_structure,type,'yz') - end - end - elseif length(tensor.xx)==md.numberofgrids | length(tensor.xx)==md.numberofgrids2d, - if ~isnan(options_structure.noiceshelf) & options_structure.noiceshelf, - pos=find(md.gridoniceshelf); - tensor.xx(pos)=NaN; - tensor.yy(pos)=NaN; - tensor.xy(pos)=NaN; - if strcmpi(md.type,'3d') - tensor.zz(pos)=NaN; - tensor.xz(pos)=NaN; - tensor.yz(pos)=NaN; - end - end - if ~isnan(options_structure.noiceshelf) & options_structure.noicesheet, - pos=find(md.gridonicesheet); - tensor.xx(pos)=NaN; - tensor.yy(pos)=NaN; - tensor.xy(pos)=NaN; - if strcmpi(md.type,'3d') - tensor.zz(pos)=NaN; - tensor.xz(pos)=NaN; - tensor.yz(pos)=NaN; - end - end - - if strcmpi(md.type,'2d'), - A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3); - - subplot(2*width,2*width,index1) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.xx ,'FaceColor','flat','EdgeColor','none'); - Apply_options_tensor(options_structure,type,'xx') - - subplot(2*width,2*width,index2) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.yy ,'FaceColor','flat','EdgeColor','none'); - Apply_options_tensor(options_structure,type,'yy') - - subplot(2*width,2*width,index3) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.xy ,'FaceColor','flat','EdgeColor','none'); - Apply_options_tensor(options_structure,type,'xy') - - else - if options_structure.layer>=1, - A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3); - subplot(3*width,3*width,index1) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.xx ,'FaceColor','flat','EdgeColor','none'); - Apply_options_tensor(options_structure,type,'xx') - - subplot(3*width,3*width,index2) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.yy ,'FaceColor','flat','EdgeColor','none'); - Apply_options_tensor(options_structure,type,'yy') - - subplot(3*width,3*width,index3) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.zz ,'FaceColor','flat','EdgeColor','none'); - Apply_options_tensor(options_structure,type,'zz') - - subplot(3*width,3*width,index4) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.xy ,'FaceColor','flat','EdgeColor','none'); - Apply_options_tensor(options_structure,type,'xy') - - subplot(3*width,3*width,index5) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.xz ,'FaceColor','flat','EdgeColor','none'); - Apply_options_tensor(options_structure,type,'xz') - - subplot(3*width,3*width,index6) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.yz ,'FaceColor','flat','EdgeColor','none'); - Apply_options_tensor(options_structure,type,'yz') - - else - A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3); D=md.elements(:,4); E=md.elements(:,5); F=md.elements(:,6); - - subplot(3*width,3*width,index1) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xx,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xx,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xx,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xx,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xx,'FaceColor','flat','EdgeColor','none'); - Apply_options_tensor(options_structure,type,'xx') - - subplot(3*width,3*width,index2) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.yy,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.yy,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.yy,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.yy,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.yy,'FaceColor','flat','EdgeColor','none'); - Apply_options_tensor(options_structure,type,'yy') - - subplot(3*width,3*width,index3) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.zz,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.zz,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.zz,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.zz,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.zz,'FaceColor','flat','EdgeColor','none'); - Apply_options_tensor(options_structure,type,'zz') - - subplot(3*width,3*width,index4) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xy,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xy,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xy,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xy,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xy,'FaceColor','flat','EdgeColor','none'); - Apply_options_tensor(options_structure,type,'xy') - - subplot(3*width,3*width,index5) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xz,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xz,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xz,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xz,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xz,'FaceColor','flat','EdgeColor','none'); - Apply_options_tensor(options_structure,type,'xz') - - subplot(3*width,3*width,index6) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.yz,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.yz,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.yz,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.yz,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.yz,'FaceColor','flat','EdgeColor','none'); - Apply_options_tensor(options_structure,type,'yz') - end - end - end +%Compute the indexes of the components plots +upperplots=fix((i-1)/width); +if upperplots==0, leftplots=i-1; else leftplots=i-width*upperplots-1; end +if strcmpi(md.type,'2d')%3 components -> 3 indexes + index1=4*width*upperplots+2*leftplots+1; + index2=index1+1; + index3=index1+width*2; +elseif strcmpi(md.type,'3d')%6 components -> 6 indexes + index1=3*3*width*upperplots+3*leftplots+1; + index2=index1+1; + index3=index1+2; + index4=index1+width*3; + index5=index4+1; + index6=index4+2; end -function Apply_options_tensor(options_structure,type,component) -%apply options +%process data and model +[x y z elements is2d]=processmesh(md,options_structure); +[tensor.xx isongrid]=processdata(md,tensor.xx,options_structure); +[tensor.yy isongrid]=processdata(md,tensor.yy,options_structure); +[tensor.xy isongrid]=processdata(md,tensor.xy,options_structure); +if strcmpi(md.type,'3d') + [tensor.xz isongrid]=processdata(md,tensor.xz,options_structure); + [tensor.yz isongrid]=processdata(md,tensor.yz,options_structure); + [tensor.zz isongrid]=processdata(md,tensor.zz,options_structure); +end + +if (strcmpi(md.type,'2d')), + subplot(2*width,2*width,index1), + plot_unit(x,y,z,elements,tensor.xx,isongrid,is2d,options_structure) + Apply_options_tensor(md,options_structure,type,'xx') + subplot(2*width,2*width,index2), + plot_unit(x,y,z,elements,tensor.yy,isongrid,is2d,options_structure) + Apply_options_tensor(md,options_structure,type,'yy') + subplot(2*width,2*width,index3), + plot_unit(x,y,z,elements,tensor.xy,isongrid,is2d,options_structure) + Apply_options_tensor(md,options_structure,type,'xy') +else + subplot(3*width,3*width,index1), + plot_unit(x,y,z,elements,tensor.xx,isongrid,is2d,options_structure) + Apply_options_tensor(md,options_structure,type,'xx') + subplot(3*width,3*width,index2), + plot_unit(x,y,z,elements,tensor.yy,isongrid,is2d,options_structure) + Apply_options_tensor(md,options_structure,type,'yy') + subplot(3*width,3*width,index3), + plot_unit(x,y,z,elements,tensor.zz,isongrid,is2d,options_structure) + Apply_options_tensor(md,options_structure,type,'zz') + subplot(3*width,3*width,index4), + plot_unit(x,y,z,elements,tensor.xy,isongrid,is2d,options_structure) + Apply_options_tensor(md,options_structure,type,'xy') + subplot(3*width,3*width,index5), + plot_unit(x,y,z,elements,tensor.xz,isongrid,is2d,options_structure) + Apply_options_tensor(md,options_structure,type,'xz') + subplot(3*width,3*width,index6), + plot_unit(x,y,z,elements,tensor.yz,isongrid,is2d,options_structure) + Apply_options_tensor(md,options_structure,type,'yz') +end +end + +function Apply_options_tensor(md,options_structure,type,component) + %apply options if isnan(options_structure.title) if ismember('_',type) %user plotet stress_tensor Index: /issm/trunk/src/m/classes/public/plot/plot_tensor_principal.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plot_tensor_principal.m (revision 26) +++ /issm/trunk/src/m/classes/public/plot/plot_tensor_principal.m (revision 27) @@ -5,211 +5,64 @@ % plot_tensor_principal(md,options_structure,width,i,tensor,type,plot_options); % -% See also: PLOTMODEL +% See also: PLOTMODEL, PLOT_UNIT, PLOT_DEALER - %Compute the indexes of the components plots - upperplots=fix((i-1)/width); - if upperplots==0, leftplots=i-1; else leftplots=i-width*upperplots-1; end - if strcmpi(md.type,'2d')%3 components -> 3 indexes - index1=4*width*upperplots+2*leftplots+1; - index2=index1+1; - index3=index1+width*2; - index4=index3+1; - newwidth=2*width; - elseif strcmpi(md.type,'3d')%6 components -> 6 indexes - index1=3*3*width*upperplots+3*leftplots+1; - index2=index1+1; - index3=index1+2; - index4=index1+width*3; - index5=index4+1; - index6=index4+2; - newwidth=3*width; - end - - %plot principal axis - type1=[type 'axis1']; - plot_tensor_principalaxis(md,options_structure,newwidth,index1,tensor,type1,plot_options); - type2=[type 'axis2']; - plot_tensor_principalaxis(md,options_structure,newwidth,index2,tensor,type2,plot_options); - if strcmpi(md.type,'3d') - type3=[type 'axis3']; - plot_tensor_principalaxis(md,options_structure,newwidth,index3,tensor,type3,plot_options); - end - - %smoothing? - if strcmpi(options_structure.smooth,'yes') & length(tensor.principalvalue1)==md.numberofelements - tensor.principalvalue1=elementstogrids(md,tensor.principalvalue1); - tensor.principalvalue2=elementstogrids(md,tensor.principalvalue2); - if strcmpi(md.type,'3d'), tensor.principalvalue3=elementstogrids(md,tensor.principalvalue3); end - end - - - %layer projection? - if ~isnan(options_structure.layer) & options_structure.layer>=1, - tensor.principalvalue1=project2d(md,tensor.principalvalue1,options_structure.layer); %project onto 2d mesh - tensor.principalvalue2=project2d(md,tensor.principalvalue2,options_structure.layer); - tensor.principalvalue3=project2d(md,tensor.principalvalue3,options_structure.layer); - - %we modify the mesh temporarily to a 2d mesh from which the 3d mesh was extruded. - md.x=md.x2d; - md.y=md.y2d; - md.z=md.z2d; - md.elements=md.elements2d; - md.elements_type=md.elements_type2d; - end - - %units - if ~isnan(options_structure.unitmultiplier), - md.x=md.x*options_structure.unitmultiplier; - md.y=md.y*options_structure.unitmultiplier; - md.z=md.z*options_structure.unitmultiplier; - end - - if length(tensor.principalvalue1)==length(md.elements), - if ~isnan(options_structure.noiceshelf) & options_structure.noiceshelf, - tensor.principalvalue1(find(md.elementoniceshelf))=NaN; - tensor.principalvalue2(find(md.elementoniceshelf))=NaN; - if strcmpi(md.type,'3d') - tensor.principalvalue3(find(md.elementoniceshelf))=NaN; - end - end - if ~isnan(options_structure.noicesheet) & options_structure.noicesheet, - tensor.principalvalue1(find(~md.elementoniceshelf))=NaN; - tensor.principalvalue2(find(~md.elementoniceshelf))=NaN; - if strcmpi(md.type,'3d') - tensor.principalvalue3(find(~md.elementoniceshelf))=NaN; - end - end - - if (strcmpi(md.type,'2d')), - A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3); - - subplot(2*width,2*width,index3) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue1,'FaceColor','flat','EdgeColor','black'); - Apply_options_tensor(options_structure,type,'principal value 1') - - subplot(2*width,2*width,index4) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue2,'FaceColor','flat','EdgeColor','black'); - Apply_options_tensor(options_structure,type,'principal value 2') - else - if options_structure.layer>=1, - A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3); - - subplot(3*width,3*width,index4) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue1,'FaceColor','flat','EdgeColor','black'); - Apply_options_tensor(options_structure,type,'principal value 1') - - subplot(3*width,3*width,index5) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue2,'FaceColor','flat','EdgeColor','black'); - Apply_options_tensor(options_structure,type,'principal value 2') - - subplot(3*width,3*width,index6) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue3,'FaceColor','flat','EdgeColor','black'); - Apply_options_tensor(options_structure,type,'principal value 3') - - else - - A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3); D=md.elements(:,4); E=md.elements(:,5); F=md.elements(:,6); - - subplot(3*width,3*width,index4) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue1,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue1,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue1,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue1,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue1,'FaceColor','flat','EdgeColor','black'); - Apply_options_tensor(options_structure,type,'principal value 1') - - subplot(3*width,3*width,index5) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue2,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue2,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue2,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue2,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue2,'FaceColor','flat','EdgeColor','black'); - Apply_options_tensor(options_structure,type,'principal value 2') - - subplot(3*width,3*width,index6) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue3,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue3,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue3,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue3,'FaceColor','flat','EdgeColor','black'); - patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue3,'FaceColor','flat','EdgeColor','black'); - Apply_options_tensor(options_structure,type,'principal value 3') - end - end - - elseif length(tensor.principalvalue1)==md.numberofgrids | length(tensor.principalvalue1)==md.numberofgrids2d, - if ~isnan(options_structure.noiceshelf) & options_structure.noiceshelf, - pos=find(md.gridoniceshelf); - tensor.principalvalue1(pos)=NaN; - tensor.principalvalue2(pos)=NaN; - if strcmpi(md.type,'3d') - tensor.principalvalue3(pos)=NaN; - end - end - if ~isnan(options_structure.noiceshelf) & options_structure.noicesheet, - pos=find(md.gridonicesheet); - tensor.principalvalue1(pos)=NaN; - tensor.principalvalue2(pos)=NaN; - if strcmpi(md.type,'3d') - tensor.principalvalue3(pos)=NaN; - end - end - - if strcmpi(md.type,'2d'), - A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3); - - subplot(2*width,2*width,index3) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.principalvalue1 ,'FaceColor','flat','EdgeColor','none'); - Apply_options_tensor(options_structure,type,'principal value 1') - - subplot(2*width,2*width,index4) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.principalvalue2 ,'FaceColor','flat','EdgeColor','none'); - Apply_options_tensor(options_structure,type,'principal value 2') - else - if options_structure.layer>=1, - A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3); - - subplot(3*width,3*width,index4) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.principalvalue1 ,'FaceColor','flat','EdgeColor','none'); - Apply_options_tensor(options_structure,type,'principal value 1') - - subplot(3*width,3*width,index5) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.principalvalue2 ,'FaceColor','flat','EdgeColor','none'); - Apply_options_tensor(options_structure,type,'principal value 2') - - subplot(3*width,3*width,index6) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.principalvalue3 ,'FaceColor','flat','EdgeColor','none'); - Apply_options_tensor(options_structure,type,'principal value 3') - else - A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3); D=md.elements(:,4); E=md.elements(:,5); F=md.elements(:,6); - - subplot(3*width,3*width,index4) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue1,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue1,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue1,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue1,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue1,'FaceColor','flat','EdgeColor','none'); - Apply_options_tensor(options_structure,type,'principal value 1') - - subplot(3*width,3*width,index5) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue2,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue2,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue2,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue2,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue2,'FaceColor','flat','EdgeColor','none'); - Apply_options_tensor(options_structure,type,'principal value 2') - - subplot(3*width,3*width,index6) - patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue3,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue3,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue3,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue3,'FaceColor','flat','EdgeColor','none'); - patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue3,'FaceColor','flat','EdgeColor','none'); - Apply_options_tensor(options_structure,type,'principal value 3') - end - end - end +%Compute the indexes of the components plots +upperplots=fix((i-1)/width); +if upperplots==0, leftplots=i-1; else leftplots=i-width*upperplots-1; end +if strcmpi(md.type,'2d')%3 components -> 3 indexes + index1=4*width*upperplots+2*leftplots+1; + index2=index1+1; + index3=index1+width*2; + index4=index3+1; + newwidth=2*width; +elseif strcmpi(md.type,'3d')%6 components -> 6 indexes + index1=3*3*width*upperplots+3*leftplots+1; + index2=index1+1; + index3=index1+2; + index4=index1+width*3; + index5=index4+1; + index6=index4+2; + newwidth=3*width; end -function Apply_options_tensor(options_structure,type,component) +%plot principal axis +type1=[type 'axis1']; +plot_tensor_principalaxis(md,options_structure,newwidth,index1,tensor,type1,plot_options); +type2=[type 'axis2']; +plot_tensor_principalaxis(md,options_structure,newwidth,index2,tensor,type2,plot_options); +if strcmpi(md.type,'3d') + type3=[type 'axis3']; + plot_tensor_principalaxis(md,options_structure,newwidth,index3,tensor,type3,plot_options); +end + +%now plot principal values +[x y z elements is2d]=processmesh(md,options_structure); +[tensor.principalvalue1 isongrid]=processdata(md,tensor.principalvalue1,options_structure); +[tensor.principalvalue2 isongrid]=processdata(md,tensor.principalvalue2,options_structure); +if strcmpi(md.type,'3d') + [tensor.principalvalue3 isongrid]=processdata(md,tensor.principalvalue3,options_structure); +end + +if (strcmpi(md.type,'2d')), + subplot(2*width,2*width,index3) + plot_unit(x,y,z,elements,tensor.principalvalue1,isongrid,is2d,options_structure) + Apply_options_tensor(md,options_structure,type,'principal value 1') + subplot(2*width,2*width,index4) + plot_unit(x,y,z,elements,tensor.principalvalue2,isongrid,is2d,options_structure) + Apply_options_tensor(md,options_structure,type,'principal value 2') +else + subplot(3*width,3*width,index4) + plot_unit(x,y,z,elements,tensor.principalvalue1,isongrid,is2d,options_structure) + Apply_options_tensor(md,options_structure,type,'principal value 1') + subplot(3*width,3*width,index5) + plot_unit(x,y,z,elements,tensor.principalvalue2,isongrid,is2d,options_structure) + Apply_options_tensor(md,options_structure,type,'principal value 2') + subplot(3*width,3*width,index6) + plot_unit(x,y,z,elements,tensor.principalvalue3,isongrid,is2d,options_structure) + Apply_options_tensor(md,options_structure,type,'principal value 3') +end +end + +function Apply_options_tensor(md,options_structure,type,component) %apply options if isnan(options_structure.title) Index: /issm/trunk/src/m/classes/public/plot/plot_tensor_principalaxis.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plot_tensor_principalaxis.m (revision 26) +++ /issm/trunk/src/m/classes/public/plot/plot_tensor_principalaxis.m (revision 27) @@ -7,87 +7,81 @@ % See also: PLOTMODEL -%plot mesh boundaries +%prepare subplot subplot(width,width,i); +%process data and model +[x y z elements is2d]=processmesh(md,options_structure); if (strcmpi(md.type,'2d')), eval(['Vx=tensor.principalaxis' type(end) '(:,1); Vy=tensor.principalaxis' type(end) '(:,2);']) + eval(['value=tensor.principalvalue' type(end) ';']); + [Vx isongrid]=processdata(md,Vx,options_structure); + [Vy isongrid]=processdata(md,Vy,options_structure); + [value isongrid]=processdata(md,value,options_structure); else eval(['Vx=tensor.principalaxis' type(end) '(:,1); Vy=tensor.principalaxis' type(end) '(:,2); Vz=tensor.principalaxis' type(end) '(:,3);']) + [Vx isongrid]=processdata(md,Vx,options_structure); + [Vy isongrid]=processdata(md,Vy,options_structure); + [Vz isongrid]=processdata(md,Vz,options_structure); + [value isongrid]=processdata(md,value,options_structure); end -%smoothing? -if strcmpi(options_structure.smooth,'yes') & length(Vx)==md.numberofelements - Vx=elementstogrids(md,Vx); - Vy=elementstogrids(md,Vy); - if (strcmpi(md.type,'3d')), - Vz=elementstogrids(md,Vz); - end -end - -%layer projection? -if ~isnan(options_structure.layer) & options_structure.layer>=1, - Vx=project2d(md,Vx,options_structure.layer); %project onto 2d mesh - Vy=project2d(md,Vy,options_structure.layer); %project onto 2d mesh - if (strcmpi(md.type,'3d')), - Vz=project2d(md,Vz,options_structure.layer); %project onto 2d mesh - end - %we modify the mesh temporarily to a 2d mesh from which the 3d mesh was extruded. - md.x=md.x2d; - md.y=md.y2d; - md.z=md.z2d; - md.elements=md.elements2d; - md.elements_type=md.elements_type2d; - md.type='2d'; +%take the center of each element if ~isongrid +if ~isongrid + x=mean(md.x(md.elements'))'; y=mean(md.y(md.elements'))'; z=mean(md.z(md.elements'))'; end -%units -if ~isnan(options_structure.unitmultiplier), - md.x=md.x*options_structure.unitmultiplier; - md.y=md.y*options_structure.unitmultiplier; - md.z=md.z*options_structure.unitmultiplier; -end -if length(Vx)==length(md.elements), - - if ~isnan(options_structure.noiceshelf) & options_structure.noiceshelf, - Vx(find(md.elementoniceshelf))=NaN; - Vy(find(md.elementoniceshelf))=NaN; - Vz(find(md.elementoniceshelf))=NaN; - end - if ~isnan(options_structure.noicesheet) & options_structure.noicesheet, - Vx(find(~md.elementoniceshelf))=NaN; - Vy(find(~md.elementoniceshelf))=NaN; - Vz(find(~md.elementoniceshelf))=NaN; +%plot quivers +if strcmpi(md.type,'2d'), + + %density + if ~isnan(options_structure.density) + x=x(1:options_structure.density:end); + y=y(1:options_structure.density:end); + Vx=Vx(1:options_structure.density:end); + Vy=Vy(1:options_structure.density:end); + value=value(1:options_structure.density:end); end - if (strcmpi(md.type,'2d')), - x=mean(md.x(md.elements'))'; y=mean(md.y(md.elements'))'; - A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3); - quiver(x,y,Vx,Vy) + %scaling: + delta=((min(x)-max(x))^2+(min(y)-max(y))^2)/numel(x); + scale=0.5/max(sqrt((Vx.^2+Vy.^2)/delta)); + Vx=scale*Vx; Vy=scale*Vy; - else - x=mean(md.x(md.elements'))'; y=mean(md.y(md.elements'))'; z=mean(md.z(md.elements'))'; - A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3); D=md.elements(:,4); E=md.elements(:,5); F=md.elements(:,6); - quiver3(x,y,z,Vx,Vy,Vz) - end -elseif length(Vx)==md.numberofgrids | length(Vx)==md.numberofgrids2d, - if ~isnan(options_structure.noiceshelf) & options_structure.noiceshelf, - pos=find(md.gridoniceshelf); - Vx(pos)=NaN; - Vy(pos)=NaN; - Vz(pos)=NaN; - end - if ~isnan(options_structure.noiceshelf) & options_structure.noicesheet, - pos=find(md.gridonicesheet); - Vx(pos)=NaN; - Vy(pos)=NaN; - Vz(pos)=NaN; + pos=find(value>=0); + q1=quiver(x(pos),y(pos),Vx(pos),Vy(pos),'Color','r','ShowArrowHead','off','AutoScale','off'); + hold on + pos=find(value<0); + q2=quiver(x(pos),y(pos),Vx(pos),Vy(pos),'Color','b','ShowArrowHead','off','AutoScale','off'); + +else + %density + if ~isnan(options_structure.density) + x=x(1:options_structure.density:end); + y=y(1:options_structure.density:end); + z=z(1:options_structure.density:end); + Vx=Vx(1:options_structure.density:end); + Vy=Vy(1:options_structure.density:end); + Vz=Vz(1:options_structure.density:end); + value=value(1:options_structure.density:end); end - if strcmpi(md.type,'2d'), - quiver(md.x,md.y,Vx,Vy) - else - quiver3(md.x,md.y,md.z,Vx,Vy,Vz) - end + %scaling: + delta=((min(x)-max(x))^2+(min(y)-max(y))^2)/numel(x); + scale=0.5/max(sqrt((Vx.^2+Vy.^2)/delta)); + Vx=scale*Vx; Vy=scale*Vy; Vz=scale*Vz; + + pos=find(value>=0); + q1=quiver3(x(pos),y(pos),z(pos),Vx(pos),Vy(pos),Vz(pos),'Color','r','ShowArrowHead','off','AutoScale','off'); + hold on + pos=find(value<0); + q2=quiver3(x(pos),y(pos),z(pos),Vx(pos),Vy(pos),Vz(pos),'Color','b','ShowArrowHead','off','AutoScale','off'); +end + +%legend +if strcmpi(type(1:6),'strain') + legend([q1 q2],'extension','compression') +elseif strcmpi(type(1:6),'stress') + legend([q1 q2],'compression','traction') end Index: /issm/trunk/src/m/classes/public/plot/plot_transient_movie.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plot_transient_movie.m (revision 26) +++ /issm/trunk/src/m/classes/public/plot/plot_transient_movie.m (revision 27) @@ -4,90 +4,40 @@ % plot_transient_movie(md,options_structure,width,i); % -% See also: PLOTMODEL +% See also: PLOTMODEL, PLOT_UNIT, PLOT_DEALER - %plot mesh boundaries + %prepare subplot subplot(width,width,i); - - %first load x,y, etc ... to speed up plot - x=md.x; - x2d=md.x2d; - y=md.y; - y2d=md.y2d; - z=md.z; - z2d=md.z2d; - elements2d=md.elements2d; - elements=md.elements; - elements_type2d=md.elements_type2d; - - %units - if ~isnan(options_structure.unitmultiplier), - md.x=md.x*options_structure.unitmultiplier; - md.y=md.y*options_structure.unitmultiplier; - md.z=md.z*options_structure.unitmultiplier; - end - - %edgecolor? - if ~isnan(options_structure.edgecolor), - edgecolor=options_structure.edgecolor; - else - edgecolor='none'; - end if strcmpi(md.type,'2d') choice=input('Which field do you want to plot? (vel/vx/vy/thickness/surface/bed)','s'); - if ~strcmp(choice,'vel') & ~strcmp(choice,'vx') & ~strcmp(choice,'vy') & ~strcmp(choice,'thickness') & ~strcmp(choice,'bed') & ~strcmp(choice,'surface') disp('plot_transient_movie error message: input not supported yet, exiting...') return end - if ~isnan(options_structure.noiceshelf) & options_structure.noiceshelf, - pos=find(md.gridoniceshelf); - data(pos)=NaN; - end - if ~isnan(options_structure.noiceshelf) & options_structure.noicesheet, - pos=find(md.gridonicesheet); - data(pos)=NaN; - end - for i=1:length(md.transient_results) - eval(['data=md.transient_results(' num2str(i) ').' num2str(choice) ';']); - titlestring=[choice ' at time ' num2str(md.transient_results(i).time) ' year']; - A=elements(:,1); B=elements(:,2); C=elements(:,3); - patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor); - apply_options_movie(options_structure,titlestring); - pause(0.5) - end - else choice=input('Which field do you want to plot? (vel/vx/vy/vz/thickness/surface/bed/temperature)','s'); - if ~strcmp(choice,'vel') & ~strcmp(choice,'vx') & ~strcmp(choice,'vy') & ~strcmp(choice,'vz') & ~strcmp(choice,'temperature') & ~strcmp(choice,'thickness') & ~strcmp(choice,'bed') & ~strcmp(choice,'surface') disp('plot_transient_movie error message: input not supported yet, exiting...') return end - if ~isnan(options_structure.noiceshelf) & options_structure.noiceshelf, - pos=find(md.gridoniceshelf); - data(pos)=NaN; - end - if ~isnan(options_structure.noiceshelf) & options_structure.noicesheet, - pos=find(md.gridonicesheet); - data(pos)=NaN; - end - for i=1:length(md.transient_results) - eval(['data=md.transient_results(' num2str(i) ').' num2str(choice) ';']); - titlestring=[choice 'at time ' num2str(md.transient_results(i).time) ' a']; - A=elements(:,1); B=elements(:,2); C=elements(:,3); D=elements(:,4); E=elements(:,5); F=elements(:,6); - patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor); - patch( 'Faces', [D E F], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor); - patch( 'Faces', [A B E D], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor); - patch( 'Faces', [B E F C ], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor); - patch( 'Faces', [C A D F ], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor); - apply_options_movie(options_structure,titlestring); - pause(0.5) - end end + %process model + [x y z elements is2d]=processmesh(md,options_structure); + + %loop over the time steps + for i=1:length(md.transient_results) + eval(['data=md.transient_results(' num2str(i) ').' num2str(choice) ';']); + + %process data + [data isongrid]=processdata(md,data,options_structure); + titlestring=[choice ' at time ' num2str(md.transient_results(i).time) ' year']; + plot_unit(x,y,z,elements,data,isongrid,is2d,options_structure) + apply_options_movie(md,options_structure,titlestring); + pause(0.5) + end end %function -function apply_options_movie(options_structure,titlestring) +function apply_options_movie(md,options_structure,titlestring) %apply options if isnan(options_structure.title) Index: /issm/trunk/src/m/classes/public/plot/plot_unit.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plot_unit.m (revision 26) +++ /issm/trunk/src/m/classes/public/plot/plot_unit.m (revision 27) @@ -1,120 +1,9 @@ -function plot_unit(md,optionstring,width,i); -%PLOT_UNIT - unit plot called by plotmodel +function plot_unit(x,y,z,elements,data,isongrid,is2d,options_structure) +%PLOT_UNIT - unit plot, display data % % Usage: -% plot_unit(md,optionstring,width,i); +% plot_unit(x,y,z,elements,data,isongrid,is2d,options_structure); % -% See also: PLOTMODEL - -%parse options and get a structure of options. -options_structure=parse_options(md,optionstring); - -%figure out if this is a special plot -data=findarg(optionstring,'data');data=data.value; - -%Figure out if this is a semi-transparent plot. -if ~isnan(options_structure.overlay), - plot_overlay(md,data,options_structure,width,i); - return; -end - -if ischar(data), - switch data, - case 'boundaries', - plot_boundaries(md,options_structure,width,i); - return; - case 'elementnumbering', - plot_elementnumbering(md,options_structure,width,i); - return; - case 'highlightelements', - plot_highlightelements(md,options_structure,width,i); - case 'segmentnumbering', - plot_segmentnumbering(md,options_structure,width,i); - return; - case 'importancefactors', - plot_importancefactors(md,options_structure,width,i); - return; - case 'elements_type', - plot_elementstype(md,options_structure,width,i); - return; - case 'gridnumbering', - plot_gridnumbering(md,options_structure,width,i); - return; - case 'highlightgrids', - plot_highlightgrids(md,options_structure,width,i); - return; - case 'basal_drag', - plot_basaldrag(md,options_structure,width,i); - return; - case 'driving_stress', - plot_drivingstress(md,options_structure,width,i); - return; - case 'mesh', - plot_mesh(md,options_structure,width,i); - return; - case 'penalties', - plot_penalties(md,options_structure,width,i); - return; - case 'quiver', - plot_quiver(md,options_structure,width,i); - return; - case 'quiver3', - plot_quiver3(md,options_structure,width,i); - return; - case 'quivervel', - plot_quivervel(md,options_structure,width,i); - return; - case 'riftvel', - plot_riftvel(md,options_structure,width,i); - return; - case 'riftrelvel', - plot_riftrelvel(md,options_structure,width,i); - return; - case 'riftpenetration', - plot_riftpenetration(md,options_structure,width,i); - return; - case 'sarpwr', - plot_sarpwr(md,options_structure,width,i) - return - case {'segmentonneumann_diag','segmentonneumann_prog'} - plot_segmentonneumann(md,options_structure,width,i,data) - return - case {'strainrate_tensor','strainrate','strainrate_principal','strainrate_principalaxis1','strainrate_principalaxis2','strainrate_principalaxis3',... - 'stress_tensor','stress','stress_principal','stress_principalaxis1','stress_principalaxis2','stress_principalaxis3',... - 'deviatoricstress_tensor','deviatoricstress','deviatoricstress_principal','deviatoricstress_principalaxis1','deviatoricstress_principalaxis2','deviatoricstress_principalaxis3'}, - plot_tensor(md,options_structure,width,i,data); - return; - case 'thermaltransient_results', - plot_thermaltransient_results(md,options_structure,width,i); - return; - case 'transient_movie', - plot_transient_movie(md,options_structure,width,i); - return; - case 'transient_results', - plot_transient_results(md,options_structure,width,i); - return; - otherwise, - if isfield(struct(md),data) - data=eval(['md.' data ';']); - else - error('plot error message: data provided not supported yet. Type plotdoc for help'); - end - end -end - -%Figure out if this is a Section plot -if ~isnan(options_structure.sectionvalue) - plot_section(md,data,options_structure,width,i); - return; -end - -%check on arguments -if (iscell(data) | isempty(data)), - error('plot error message: data provided is empty'); -end - -%process data and model -[x y z elements is2d]=processmesh(md,options_structure); -[data isongrid]=processdata(md,data,options_structure); +% See also: PLOTMODEL, PLOT_DEALER %edgecolor? @@ -124,7 +13,4 @@ edgecolor='none'; end - -%standard plot: -subplot(width,width,i); %element data @@ -145,24 +31,17 @@ if is2d A=elements(:,1); B=elements(:,2); C=elements(:,3); - patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor); + patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', data(:),'FaceColor','interp','EdgeColor',edgecolor); else if options_structure.layer>=1, A=elements(:,1); B=elements(:,2); C=elements(:,3); - patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor); + patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', data(:),'FaceColor','interp','EdgeColor',edgecolor); else A=elements(:,1); B=elements(:,2); C=elements(:,3); D=elements(:,4); E=elements(:,5); F=elements(:,6); - patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor); - patch( 'Faces', [D E F], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor); - patch( 'Faces', [A B E D], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor); - patch( 'Faces', [B E F C ], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor); - patch( 'Faces', [C A D F ], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor); + patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', data(:),'FaceColor','interp','EdgeColor',edgecolor); + patch( 'Faces', [D E F], 'Vertices', [x y z],'FaceVertexCData', data(:),'FaceColor','interp','EdgeColor',edgecolor); + patch( 'Faces', [A B E D], 'Vertices', [x y z],'FaceVertexCData', data(:),'FaceColor','interp','EdgeColor',edgecolor); + patch( 'Faces', [B E F C ], 'Vertices', [x y z],'FaceVertexCData', data(:),'FaceColor','interp','EdgeColor',edgecolor); + patch( 'Faces', [C A D F ], 'Vertices', [x y z],'FaceVertexCData', data(:),'FaceColor','interp','EdgeColor',edgecolor); end end end - -%apply all options -if isnan(options_structure.shading) & isnan(options_structure.edgecolor) & size(data,1)==md.numberofgrids, - options_structure.shading='interp'; -end - -applyoptions(md,data,options_structure); Index: /issm/trunk/src/m/classes/public/plot/plotdoc.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plotdoc.m (revision 26) +++ /issm/trunk/src/m/classes/public/plot/plotdoc.m (revision 27) @@ -30,4 +30,5 @@ disp(' - ''riftpenetration'': penetration levels for a fault'); disp(' - ''basal_drag'': plot the basal drag on the bed (in kPa)'); +disp(' - ''basal_dragx'' or ''basal_dragy'' : plot a component of the basal drag on the bed (in kPa)'); disp(' - ''driving_stress'': plot the driving stress (in kPa)'); disp(' - ''strainrate_tensor'': plot the components of the strainrate tensor (exx,eyy,ezz,exy,exz,eyz) if computed'); @@ -53,5 +54,7 @@ disp(' ''contourticks'': ''on'' or ''off'' to display the ticks of the contours'); disp(' ''contouronly'': ''on'' or ''off'' to display the contours on a white background'); -disp(' ''contourcolors'': ticks and contour color'); +disp(' ''contourcolor'': ticks and contour color'); +disp(' ''density'': density of quivers (one arrow every N nodes, N integer)'); +disp(' ''streamlines'': N (number of stream lines) or {[x1 y1],...} (coordinates of seed points) add streanlines on current figure'); disp(' ''wrapping'': repeat ''n'' times the colormap (''n'' must be an integer)'); disp(' ''edgecolor'': same as standard matlab option EdgeColor (color name: ''black'' or RGB array: [0.5 0.2 0.8])'); @@ -61,5 +64,5 @@ disp(' ''resolution'': resolution used by section value (array of type [horizontal_resolution vertical_resolution])'); disp(' horizontal_resolution must be in meter, and vertical_resolution a number of layers'); -disp(' ''showsection'': show section used by ''sectionvalue'' (string ''yes'')'); +disp(' ''showsection'': show section used by ''sectionvalue'' (string ''on'' or a number of labels)'); disp(' ''sectionvalue'': give the value of data on a profile given by an Argus file (string ''Argusfile_name.exp'')'); disp(' ''smooth'': smooth element data (string ''yes'' or integer)'); Index: /issm/trunk/src/m/classes/public/plot/plotmodel.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plotmodel.m (revision 26) +++ /issm/trunk/src/m/classes/public/plot/plotmodel.m (revision 27) @@ -1,4 +1,5 @@ function plotmodel(md,varargin) %At command prompt, type plotdoc for help on documentation. + global ISSM_DIR if isempty(ISSM_DIR), @@ -28,14 +29,14 @@ %recover options for all subplots. options=recover_plot_options(md,varargin,numberofplots); + + %Create figure + figure(figurenumber),clf; - %Create figure - figure(figurenumber); clf; - %Go through all data plottable for i=1:numberofplots, - + %call unit plot - plot_unit(md,options{i},subplotwidth,i); - + plot_dealer(md,options{i},subplotwidth,i); + end else Index: /issm/trunk/src/m/classes/public/plot/processdata.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/processdata.m (revision 26) +++ /issm/trunk/src/m/classes/public/plot/processdata.m (revision 27) @@ -7,10 +7,15 @@ % See also: PLOTMODEL, PROCESSMESH -%some checks -if length(data)~=md.numberofgrids & length(data)~=md.numberofelements & length(data)~=md.numberofgrids*6 +%check format +if (iscell(data) | isempty(data)), + error('plotmodel error message: data provided is empty'); +end + +%check length +if length(data)~=md.numberofgrids & length(data)~=md.numberofelements & length(data)~=md.numberofgrids*6 & (strcmpi(md.type,'2d') | length(data)~=md.numberofgrids2d) error('plotmodel error message: data not supproted yet') end -%6*grid data +%treat the case length(data)=6*grids if length(data)==6*md.numberofgrids %keep the only norm of data @@ -18,4 +23,10 @@ data2=data(2:6:md.numberofgrids*6); data=sqrt(data1.^2+data2.^2); + %---> go to grid data +end + +%treat the case length(data)=grids2d +if (strcmpi(md.type,'3d') & length(data)==md.numberofgrids2d), + data=project3d(md,data,'node'); %---> go to grid data end @@ -60,4 +71,2 @@ data=project2d(md,data,options_structure.layer); %project onto 2d mesh end - - Index: /issm/trunk/src/m/classes/public/solve.m =================================================================== --- /issm/trunk/src/m/classes/public/solve.m (revision 26) +++ /issm/trunk/src/m/classes/public/solve.m (revision 27) @@ -14,5 +14,5 @@ %some checks on list of arguments -solutions={'mesh2grid','mesh','thermaltransient','thermalsteady','qmu','diagnostic','diagnostic_horiz','prognostic','transient','paramters','control'}; +solutions={'mesh2grid','mesh','thermaltransient','thermalsteady','qmu','diagnostic','diagnostic_horiz','prognostic','transient','parameters','control'}; found=0; for i=1:length(solutions), @@ -24,5 +24,5 @@ if found==0, - error(['solve error messae: solution type ' solutiontype ' not supported yet!']); + error(['solve error message: solution type ' solutiontype ' not supported yet!']); end Index: /issm/trunk/src/m/classes/public/thicknessevolution.m =================================================================== --- /issm/trunk/src/m/classes/public/thicknessevolution.m (revision 26) +++ /issm/trunk/src/m/classes/public/thicknessevolution.m (revision 27) @@ -4,5 +4,5 @@ % This routine compute the new thickness of a model after a time step % according to the following formula: -% dh/dt=div(Hu)+Ms-Mb +% dh/dt=div(Hu) % % Usage: @@ -13,19 +13,30 @@ end +%load some variables (it is much faster if the variab;es are loaded from md once for all) +numberofelements=md.numberofelements; +H=md.thickness; +vx=md.vx; +vy=md.vy; +index=md.elements; +x=md.x; y=md.y; + +%initialization alpha=zeros(md.numberofelements,3); beta=zeros(md.numberofelements,3); +gradx=zeros(md.numberofgrids,1); +grady=zeros(md.numberofgrids,1); -for n=1:md.numberofelements - X=inv([md.x(md.elements(n,:)) md.y(md.elements(n,:)) ones(3,1)]); - alpha(n,:)=X(1,:); - beta(n,:)=X(2,:); -end +%build some usefull variables +line=index(:); +summation=1/3*ones(3,1); +linesize=3*numberofelements; +x1=x(index(:,1)); x2=x(index(:,2)); x3=x(index(:,3)); y1=y(index(:,1)); y2=y(index(:,2)); y3=y(index(:,3)); -summation=[1;1;1]; -dHu=(md.vx(md.elements).*md.thickness(md.elements).*alpha)*summation; -dHv=(md.vy(md.elements).*md.thickness(md.elements).*beta)*summation; +%compute nodal functions coefficients N(x,y)=alpha x + beta y + gamma +invdet=1./(x1.*(y2-y3)-x2.*(y1-y3)+x3.*(y1-y2)); +alpha=[invdet.*(y2-y3) invdet.*(y3-y1) invdet.*(y1-y2)]; +beta=[invdet.*(x3-x2) invdet.*(x1-x3) invdet.*(x2-x1)]; -dHu=averaging(md,dHu,0); -dHv=averaging(md,dHv,0); - -dhdt=dHu+dHv+md.accumulation-md.melting; +%compute dhdt=div(Hu) +Hvx=H.*vx; Hvy=H.*vy; +dhdt=sum(Hvx(index).*alpha,2)+sum(Hvy(index).*beta,2);