source: issm/trunk/src/m/solutions/jpl/diagnostic.m@ 3203

function md=diagnostic(md);
%DIAGNOSTIC - compute the velocity field of a model
%
%   Usage:
%      md=diagnostic(md)
%
        %timing
        t1=clock;

        %Build all models requested for diagnostic simulation
        models.analysis_type=DiagnosticAnalysisEnum; %needed for processresults

        displaystring(md.verbose,'%s',['reading diagnostic horiz model data']);
        md.analysis_type=DiagnosticAnalysisEnum; md.sub_analysis_type=HorizAnalysisEnum; models.dh=CreateFemModel(md);

        displaystring(md.verbose,'\n%s',['reading diagnostic vert model data']);
        md.analysis_type=DiagnosticAnalysisEnum; md.sub_analysis_type=VertAnalysisEnum; models.dv=CreateFemModel(md);
        
        displaystring(md.verbose,'\n%s',['reading diagnostic stokes model data']);
        md.analysis_type=DiagnosticAnalysisEnum; md.sub_analysis_type=StokesAnalysisEnum; models.ds=CreateFemModel(md);
        
        displaystring(md.verbose,'\n%s',['reading diagnostic hutter model data']);
        md.analysis_type=DiagnosticAnalysisEnum; md.sub_analysis_type=HutterAnalysisEnum; models.dhu=CreateFemModel(md);
        
        displaystring(md.verbose,'\n%s',['reading surface and bed slope computation model data']);
        md.analysis_type=SlopecomputeAnalysisEnum; md.sub_analysis_type=NoneAnalysisEnum; models.sl=CreateFemModel(md);
        
        % figure out number of dof: just for information purposes.
        md.dof=modelsize(models);

        %initialize inputs
        inputs=inputlist;
        inputs=add(inputs,'velocity',models.dh.parameters.u_g,'doublevec',3,models.dh.parameters.numberofnodes);
        if md.control_analysis,
                inputs=add(inputs,'velocity_obs',models.dh.parameters.u_g_obs,'doublevec',2,models.dh.parameters.numberofnodes);
                inputs=add(inputs,'weights',models.dh.parameters.weights,'doublevec',1,models.dh.parameters.numberofnodes);
        end

        %compute solution
        if ~models.dh.parameters.qmu_analysis,
                if md.control_analysis,
                        %launch core of control solution.
                        results=control_core(models,inputs);

                        %process results
                        if ~isstruct(md.results), md.results=struct(); end
                        md.results.diagnostic=processresults(models,results);
                else,
                        %launch core of diagnostic solution.
                        results=diagnostic_core(models,inputs);

                        %process results
                        if ~isstruct(md.results), md.results=struct(); end
                        md.results.diagnostic=processresults(models,results);
                end
        else
                %launch dakota driver for diagnostic core solution
                Qmu(models,inputs,models.dh.parameters);
        end

        %stop timing
        t2=clock;
        displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']);