Changeset 27


Ignore:
Timestamp:
04/24/09 10:14:10 (16 years ago)
Author:
Mathieu Morlighem
Message:

updated classes form ice1

Location:
issm/trunk/src/m/classes
Files:
14 added
1 deleted
42 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/m/classes/@continuous_design/continuous_design.m

    r1 r27  
    77    properties
    88        descriptor='';
    9         initpt    = NaN;
    10         lower     = NaN;
    11         upper     = NaN;
     9        initpt    = 0.;
     10        lower     =-Inf;
     11        upper     = Inf;
     12        scale_type='none';
     13        scale     = 1.;
    1214    end
    1315   
    … …  
    3537                otherwise
    3638                    cdv.descriptor=varargin{1};
    37                     cdv.initpt    =varargin{2};
    38                     if (nargin >= 3)
    39                         cdv.lower     =varargin{3};
    40                         if (nargin >= 4)
    41                             cdv.upper     =varargin{4};
    42                             if (nargin > 4)
    43                                 warning('continuous_design:extra_arg',...
    44                                     'Extra arguments for object of class ''%s''.',...
    45                                     class(cdv));
     39
     40                    if (nargin >= 2)
     41                        cdv.initpt    =varargin{2};
     42                        if (nargin >= 3)
     43                            cdv.lower     =varargin{3};
     44                            if (nargin >= 4)
     45                                cdv.upper     =varargin{4};
     46                                if (nargin >= 5)
     47                                    cdv.scale_type=varargin{5};
     48                                    if (nargin >= 6)
     49                                        cdv.scale     =varargin{6};
     50                                        if (nargin > 6)
     51                                            warning('continuous_design:extra_arg',...
     52                                                'Extra arguments for object of class ''%s''.',...
     53                                                class(cdv));
     54                                        end
     55                                    end
     56                                end
    4657                            end
    4758                        end
    … …  
    5061
    5162        end
    52         function [desc]  =dvar_desc(cdv)
     63        function [desc]  =prop_desc(cdv)
    5364            desc=cell(size(cdv));
    5465            for i=1:numel(cdv)
    … …  
    5768            desc=allempty(desc);
    5869        end
    59         function [initpt]=dvar_initpt(cdv)
     70        function [initpt]=prop_initpt(cdv)
    6071            initpt=zeros(size(cdv));
    6172            for i=1:numel(cdv)
    6273                initpt(i)=cdv(i).initpt;
    6374            end
     75            initpt=allequal(initpt,0.);
    6476        end
    65         function [lower] =dvar_lower(cdv)
     77        function [lower] =prop_lower(cdv)
    6678            lower=zeros(size(cdv));
    6779            for i=1:numel(cdv)
    6880                lower(i)=cdv(i).lower;
    6981            end
    70             lower=allnan(lower);
     82            lower=allequal(lower,-Inf);
    7183        end
    72         function [upper] =dvar_upper(cdv)
     84        function [upper] =prop_upper(cdv)
    7385            upper=zeros(size(cdv));
    7486            for i=1:numel(cdv)
    7587                upper(i)=cdv(i).upper;
    7688            end
    77             upper=allnan(upper);
     89            upper=allequal(upper, Inf);
    7890        end
    79         function [mean]  =dvar_mean(cdv)
     91        function [mean]  =prop_mean(cdv)
    8092            mean=[];
    8193        end
    82         function [stddev]=dvar_stddev(cdv)
     94        function [stddev]=prop_stddev(cdv)
    8395            stddev=[];
    8496        end
    85         function [initst]=dvar_initst(cdv)
     97        function [initst]=prop_initst(cdv)
    8698            initst=[];
     99        end
     100        function [stype] =prop_stype(cdv)
     101            stype=cell(size(cdv));
     102            for i=1:numel(cdv)
     103                stype(i)=cellstr(cdv(i).scale_type);
     104            end
     105            stype=allequal(stype,'none');
     106        end
     107        function [scale] =prop_scale(cdv)
     108            scale=zeros(size(cdv));
     109            for i=1:numel(cdv)
     110                scale(i)=cdv(i).scale;
     111            end
     112            scale=allequal(scale,1.);
    87113        end
    88114    end

  • issm/trunk/src/m/classes/@continuous_design/display.m

    r1 r27  
    2020    disp(sprintf('        initpt: %g'      ,cdv(i).initpt));
    2121    disp(sprintf('         lower: %g'      ,cdv(i).lower));
    22     disp(sprintf('         upper: %g\n'    ,cdv(i).upper));
     22    disp(sprintf('         upper: %g'      ,cdv(i).upper));
     23    disp(sprintf('    scale_type: ''%s'''  ,cdv(i).scale_type));
     24    disp(sprintf('         scale: %g'      ,cdv(i).scale));
    2325end
    2426

  • issm/trunk/src/m/classes/@continuous_state/continuous_state.m

    r1 r27  
    77    properties
    88        descriptor='';
    9         initst    = NaN;
    10         lower     = NaN;
    11         upper     = NaN;
     9        initst    = 0.;
     10        lower     =-Inf;
     11        upper     = Inf;
    1212    end
    1313   
    … …  
    3535                otherwise
    3636                    csv.descriptor=varargin{1};
    37                     csv.initst    =varargin{2};
    38                     if (nargin >= 3)
    39                         csv.lower     =varargin{3};
    40                         if (nargin >= 4)
    41                             csv.upper     =varargin{4};
    42                             if (nargin > 4)
    43                                 warning('continuous_state:extra_arg',...
    44                                     'Extra arguments for object of class ''%s''.',...
    45                                     class(csv));
     37
     38                    if (nargin >= 2)
     39                        csv.initst    =varargin{2};
     40                        if (nargin >= 3)
     41                            csv.lower     =varargin{3};
     42                            if (nargin >= 4)
     43                                csv.upper     =varargin{4};
     44                                if (nargin > 4)
     45                                    warning('continuous_state:extra_arg',...
     46                                        'Extra arguments for object of class ''%s''.',...
     47                                        class(csv));
     48                                end
    4649                            end
    4750                        end
    … …  
    5053
    5154        end
    52         function [desc]  =dvar_desc(csv)
     55        function [desc]  =prop_desc(csv)
    5356            desc=cell(size(csv));
    5457            for i=1:numel(csv)
    … …  
    5760            desc=allempty(desc);
    5861        end
    59         function [initpt]=dvar_initpt(csv)
     62        function [initpt]=prop_initpt(csv)
    6063            initpt=[];
    6164        end
    62         function [lower] =dvar_lower(csv)
     65        function [lower] =prop_lower(csv)
    6366            lower=zeros(size(csv));
    6467            for i=1:numel(csv)
    6568                lower(i)=csv(i).lower;
    6669            end
    67             lower=allnan(lower);
     70            lower=allequal(lower,-Inf);
    6871        end
    69         function [upper] =dvar_upper(csv)
     72        function [upper] =prop_upper(csv)
    7073            upper=zeros(size(csv));
    7174            for i=1:numel(csv)
    7275                upper(i)=csv(i).upper;
    7376            end
    74             upper=allnan(upper);
     77            upper=allequal(upper, Inf);
    7578        end
    76         function [mean]  =dvar_mean(csv)
     79        function [mean]  =prop_mean(csv)
    7780            mean=[];
    7881        end
    79         function [stddev]=dvar_stddev(csv)
     82        function [stddev]=prop_stddev(csv)
    8083            stddev=[];
    8184        end
    82         function [initst]=dvar_initst(csv)
     85        function [initst]=prop_initst(csv)
    8386            initst=zeros(size(csv));
    8487            for i=1:numel(csv)
    8588                initst(i)=csv(i).initst;
    8689            end
     90            initst=allequal(initst,0.);
     91        end
     92        function [stype] =prop_stype(csv)
     93            stype={};
     94        end
     95        function [scale] =prop_scale(csv)
     96            scale=[];
    8797        end
    8898    end

  • issm/trunk/src/m/classes/@least_squares_term/least_squares_term.m

    r1 r27  
    4848
    4949        end
    50         function [desc]  =dresp_desc(lst)
     50        function [desc]  =prop_desc(lst)
    5151            desc=cell(size(lst));
    5252            for i=1:numel(lst)
    5353                desc(i)=cellstr(lst(i).descriptor);
    5454            end
     55            desc=allempty(desc);
    5556        end
    56         function [stype ]=dresp_stype(lst)
     57        function [stype] =prop_stype(lst)
    5758            stype=cell(size(lst));
    5859            for i=1:numel(lst)
    5960                stype(i)=cellstr(lst(i).scale_type);
    6061            end
     62            stype=allequal(stype,'none');
    6163        end
    62         function [scale] =dresp_scale(lst)
     64        function [scale] =prop_scale(lst)
    6365            scale=zeros(size(lst));
    6466            for i=1:numel(lst)
    6567                scale(i)=lst(i).scale;
    6668            end
     69            scale=allequal(scale,1.);
    6770        end
    68         function [weight]=dresp_weight(lst)
     71        function [weight]=prop_weight(lst)
    6972            weight=zeros(size(lst));
    7073            for i=1:numel(lst)
    7174                weight(i)=lst(i).weight;
    7275            end
     76            weight=allequal(weight,1.);
    7377        end
    74         function [lower] =dresp_lower(lst)
     78        function [lower] =prop_lower(lst)
    7579            lower=[];
    7680        end
    77         function [upper] =dresp_upper(lst)
     81        function [upper] =prop_upper(lst)
    7882            upper=[];
    7983        end
    80         function [target]=dresp_target(lst)
     84        function [target]=prop_target(lst)
    8185            target=[];
    8286        end

  • issm/trunk/src/m/classes/@model/model.m

    r1 r27  
    8888       
    8989        %Materials parameters
    90         md.rho_ice=0;
    91         md.rho_water=0;
     90        md.rho_ice=917;
     91        md.rho_water=1023;
    9292        md.heatcapacity=2093;
    9393        md.latentheat=3.34*10^5; %(J/kg);
    … …  
    9999       
    100100        %Physical parameters
    101         md.g=0;
     101        md.g=9.81;
    102102        md.yts=365*24*3600;
    103103        md.drag=NaN;
    … …  
    148148
    149149        %Statics parameters
    150         md.eps_rel=0;
    151         md.eps_abs=0;
     150        md.eps_rel=0.01;
     151        md.eps_abs=10;
    152152        md.acceleration=0;
    153         md.sparsity=0;
     153        md.sparsity=0.001;
    154154        md.connectivity=10;
    155155        md.lowmem=0;
    … …  
    172172
    173173        %Control
    174         md.control_type='drag';
     174        md.control_type={'drag'};
    175175        md.cont_vx=NaN;
    176176        md.cont_vy=NaN;
    … …  
    181181        md.nsteps=0;
    182182        md.maxiter=[];
    183         md.tolx=0;
     183        md.tolx=10^-4;
    184184        md.optscal=[];
    185185        md.mincontrolconstraint=0;
    … …  
    188188        md.epsvel=eps;
    189189        md.meanvel=0;
    190        
     190
    191191        %Output parameters
    192192        md.parameteroutput={};
    … …  
    199199        md.strainrate=NaN;
    200200        md.plot=1;
    201        
     201
    202202        %debugging
    203203        md.debug=1;
    … …  
    247247        %Cielo parameters
    248248        md.solverstring=' -mat_type aijmumps -ksp_type preonly -pc_type lu -mat_mumps_icntl_14 40 ';
    249        
    250         %Needed for running in parallel
     249
     250        %Trash, still to organize
    251251        md.analysis_type='';
    252252
    … …  
    260260
    261261        %qmu
    262         md.variables=struct([]);
    263         md.responses=struct([]);
     262        md.variables =struct();
     263        md.responses =struct();
     264        md.qmu_method=struct();
     265        md.qmu_params=struct();
    264266        md.dakotaresults=NaN;
    265         md.dakotain =NaN;
    266         md.dakotaout=NaN;
    267         md.dakotadat=NaN;
    268 %       md.method               ='sampling';
    269         md.seed                 =1234;
    270         md.samples              =10;
    271         md.sample_type          ='lhs';
    272         md.method               ='local reliability';
    273         md.interval_type        ='forward';
    274         md.fd_gradient_step_size= 0.01;
    275         md.evaluation_concurrency= 2;
    276         md.npart                 = 1;
    277 
    278         md.analysis_driver    ='';
    279         md.analysis_components='';
     267        md.dakotain     =NaN;
     268        md.dakotaout    =NaN;
     269        md.dakotadat    =NaN;
     270
     271        md.npart   =1;
    280272
    281273        %Solver options

  • issm/trunk/src/m/classes/@nonlinear_equality_constraint/nonlinear_equality_constraint.m

    r1 r27  
    77    properties
    88        descriptor='';
    9         target    = NaN;
     9        target    = 0.;
    1010        scale_type='none';
    1111        scale     = 1.;
    … …  
    5252
    5353        end
    54         function [desc]  =dresp_desc(nec)
     54        function [desc]  =prop_desc(nec)
    5555            desc=cell(size(nec));
    5656            for i=1:numel(nec)
    5757                desc(i)=cellstr(nec(i).descriptor);
    5858            end
     59            desc=allempty(desc);
    5960        end
    60         function [stype ]=dresp_stype(nec)
     61        function [stype] =prop_stype(nec)
    6162            stype=cell(size(nec));
    6263            for i=1:numel(nec)
    6364                stype(i)=cellstr(nec(i).scale_type);
    6465            end
     66            stype=allequal(stype,'none');
    6567        end
    66         function [scale] =dresp_scale(nec)
     68        function [scale] =prop_scale(nec)
    6769            scale=zeros(size(nec));
    6870            for i=1:numel(nec)
    6971                scale(i)=nec(i).scale;
    7072            end
     73            scale=allequal(scale,1.);
    7174        end
    72         function [weight]=dresp_weight(nec)
     75        function [weight]=prop_weight(nec)
    7376            weight=[];
    7477        end
    75         function [lower] =dresp_lower(nec)
     78        function [lower] =prop_lower(nec)
    7679            lower=[];
    7780        end
    78         function [upper] =dresp_upper(nec)
     81        function [upper] =prop_upper(nec)
    7982            upper=[];
    8083        end
    81         function [target]=dresp_target(nec)
     84        function [target]=prop_target(nec)
    8285            target=zeros(size(nec));
    8386            for i=1:numel(nec)
    8487                target(i)=nec(i).target;
    8588            end
     89            target=allequal(target,0.);
    8690        end
    8791    end

  • issm/trunk/src/m/classes/@nonlinear_inequality_constraint/nonlinear_inequality_constraint.m

    r1 r27  
    77    properties
    88        descriptor='';
    9         lower     = NaN;
    10         upper     = NaN;
     9        lower     =-Inf;
     10        upper     = 0.;
    1111        scale_type='none';
    1212        scale     = 1.;
    … …  
    5454
    5555        end
    56         function [desc]  =dresp_desc(nic)
     56        function [desc]  =prop_desc(nic)
    5757            desc=cell(size(nic));
    5858            for i=1:numel(nic)
    5959                desc(i)=cellstr(nic(i).descriptor);
    6060            end
     61            desc=allempty(desc);
    6162        end
    62         function [stype ]=dresp_stype(nic)
     63        function [stype] =prop_stype(nic)
    6364            stype=cell(size(nic));
    6465            for i=1:numel(nic)
    6566                stype(i)=cellstr(nic(i).scale_type);
    6667            end
     68            stype=allequal(stype,'none');
    6769        end
    68         function [scale] =dresp_scale(nic)
     70        function [scale] =prop_scale(nic)
    6971            scale=zeros(size(nic));
    7072            for i=1:numel(nic)
    7173                scale(i)=nic(i).scale;
    7274            end
     75            scale=allequal(scale,1.);
    7376        end
    74         function [weight]=dresp_weight(nic)
     77        function [weight]=prop_weight(nic)
    7578            weight=[];
    7679        end
    77         function [lower] =dresp_lower(nic)
     80        function [lower] =prop_lower(nic)
    7881            lower=zeros(size(nic));
    7982            for i=1:numel(nic)
    8083                lower(i)=nic(i).lower;
    8184            end
     85            lower=allequal(lower,-Inf);
    8286        end
    83         function [upper] =dresp_upper(nic)
     87        function [upper] =prop_upper(nic)
    8488            upper=zeros(size(nic));
    8589            for i=1:numel(nic)
    8690                upper(i)=nic(i).upper;
    8791            end
     92            upper=allequal(upper,0.);
    8893        end
    89         function [target]=dresp_target(nic)
     94        function [target]=prop_target(nic)
    9095            target=[];
    9196        end

  • issm/trunk/src/m/classes/@normal_uncertain/normal_uncertain.m

    r1 r27  
    99        mean      = NaN;
    1010        stddev    = NaN;
    11         lower     = NaN;
    12         upper     = NaN;
     11        lower     =-Inf;
     12        upper     = Inf;
    1313    end
    1414   
    … …  
    4444                    nuv.mean      =varargin{2};
    4545                    nuv.stddev    =varargin{3};
     46
    4647                    if (nargin >= 4)
    4748                        nuv.lower     =varargin{4};
    … …  
    5859
    5960        end
    60         function [desc]  =dvar_desc(nuv)
     61        function [desc]  =prop_desc(nuv)
    6162            desc=cell(size(nuv));
    6263            for i=1:numel(nuv)
    … …  
    6566            desc=allempty(desc);
    6667        end
    67         function [initpt]=dvar_initpt(nuv)
     68        function [initpt]=prop_initpt(nuv)
    6869            initpt=[];
    6970        end
    70         function [lower] =dvar_lower(nuv)
     71        function [lower] =prop_lower(nuv)
    7172            lower=zeros(size(nuv));
    7273            for i=1:numel(nuv)
    7374                lower(i)=nuv(i).lower;
    7475            end
    75             lower=allnan(lower);
     76            lower=allequal(lower,-Inf);
    7677        end
    77         function [upper] =dvar_upper(nuv)
     78        function [upper] =prop_upper(nuv)
    7879            upper=zeros(size(nuv));
    7980            for i=1:numel(nuv)
    8081                upper(i)=nuv(i).upper;
    8182            end
    82             upper=allnan(upper);
     83            upper=allequal(upper, Inf);
    8384        end
    84         function [mean]  =dvar_mean(nuv)
     85        function [mean]  =prop_mean(nuv)
    8586            mean=zeros(size(nuv));
    8687            for i=1:numel(nuv)
    … …  
    8889            end
    8990        end
    90         function [stddev]=dvar_stddev(nuv)
     91        function [stddev]=prop_stddev(nuv)
    9192            stddev=zeros(size(nuv));
    9293            for i=1:numel(nuv)
    … …  
    9495            end
    9596        end
    96         function [initst]=dvar_initst(nuv)
     97        function [initst]=prop_initst(nuv)
    9798            initst=[];
     99        end
     100        function [stype] =prop_stype(nuv)
     101            stype={};
     102        end
     103        function [scale] =prop_scale(nuv)
     104            scale=[];
    98105        end
    99106    end

  • issm/trunk/src/m/classes/@objective_function/objective_function.m

    r1 r27  
    4848
    4949        end
    50         function [desc]  =dresp_desc(of)
     50        function [desc]  =prop_desc(of)
    5151            desc=cell(size(of));
    5252            for i=1:numel(of)
    5353                desc(i)=cellstr(of(i).descriptor);
    5454            end
     55            desc=allempty(desc);
    5556        end
    56         function [stype ]=dresp_stype(of)
     57        function [stype] =prop_stype(of)
    5758            stype=cell(size(of));
    5859            for i=1:numel(of)
    5960                stype(i)=cellstr(of(i).scale_type);
    6061            end
     62            stype=allequal(stype,'none');
    6163        end
    62         function [scale] =dresp_scale(of)
     64        function [scale] =prop_scale(of)
    6365            scale=zeros(size(of));
    6466            for i=1:numel(of)
    6567                scale(i)=of(i).scale;
    6668            end
     69            scale=allequal(scale,1.);
    6770        end
    68         function [weight]=dresp_weight(of)
     71        function [weight]=prop_weight(of)
    6972            weight=zeros(size(of));
    7073            for i=1:numel(of)
    7174                weight(i)=of(i).weight;
    7275            end
     76            weight=allequal(weight,1.);
    7377        end
    74         function [lower] =dresp_lower(of)
     78        function [lower] =prop_lower(of)
    7579            lower=[];
    7680        end
    77         function [upper] =dresp_upper(of)
     81        function [upper] =prop_upper(of)
    7882            upper=[];
    7983        end
    80         function [target]=dresp_target(of)
     84        function [target]=prop_target(of)
    8185            target=[];
    8286        end

  • issm/trunk/src/m/classes/@response_function/response_function.m

    r1 r27  
    5252
    5353        end
    54         function [desc]  =dresp_desc(rf)
     54        function [desc]  =prop_desc(rf)
    5555            desc=cell(size(rf));
    5656            for i=1:numel(rf)
    5757                desc(i)=cellstr(rf(i).descriptor);
    5858            end
     59            desc=allempty(desc);
    5960        end
    60         function [stype ]=dresp_stype(rf)
     61        function [stype] =prop_stype(rf)
    6162            stype={};
    6263        end
    63         function [scale] =dresp_scale(rf)
     64        function [scale] =prop_scale(rf)
    6465            scale=[];
    6566        end
    66         function [weight]=dresp_weight(rf)
     67        function [weight]=prop_weight(rf)
    6768            weight=[];
    6869        end
    69         function [lower] =dresp_lower(rf)
     70        function [lower] =prop_lower(rf)
    7071            lower=[];
    7172        end
    72         function [upper] =dresp_upper(rf)
     73        function [upper] =prop_upper(rf)
    7374            upper=[];
    7475        end
    75         function [target]=dresp_target(rf)
     76        function [target]=prop_target(rf)
    7677            target=[];
    7778        end
    78         function [respl,probl,rell,grell]=dresp_levels(rf)
     79        function [respl,probl,rell,grell]=prop_levels(rf)
    7980            respl=cell(size(rf));
    8081            probl=cell(size(rf));

  • issm/trunk/src/m/classes/public/SectionValues.m

    r1 r27  
    8383        %Get bed and surface for each 2d point, offset to make sure that it is inside the glacier system
    8484        offset=10^-10;
    85         bed=griddata_mesh_to_mesh_3d(md.elements,md.x,md.y,md.z,md.bed,X,Y,Z,'node')+offset;
    86         surface=griddata_mesh_to_mesh_3d(md.elements,md.x,md.y,md.z,md.surface,X,Y,Z,'node')-offset;
     85        bed=griddata_mesh_to_mesh(md.elements2d,md.x2d,md.y2d,project2d(md,md.bed,1),X,Y,'node')+offset;
     86        surface=griddata_mesh_to_mesh(md.elements2d,md.x2d,md.y2d,project2d(md,md.surface,1),X,Y,'node')+offset;
    8787
    8888        %Some useful parameters

  • issm/trunk/src/m/classes/public/ThicknessCorrection.m

    r1 r27  
    2020
    2121in=ArgusContourToMesh(md.elements,md.x,md.y,expread(filename,1),'node',1);
    22 pos_shelf=find(in & md.gridoniceshelf);
    23 pos_sheet=find(in & ~md.gridoniceshelf);
     22pos_shelf=find(in & ~md.gridonicesheet);
     23pos_sheet=find(in & md.gridonicesheet);
    2424for i=1:length(pos_shelf)
    2525        %search the grid on ice sheet the closest to i

  • issm/trunk/src/m/classes/public/displaycontrol.m

    r1 r27  
    1010
    1111disp(sprintf('      ''%s''','control'));
    12 disp(sprintf('         control_type: %s (control type, ex: ''drag'', or ''B'')',md.control_type));
     12%control type
     13control_string='';
     14for i=1:length(md.control_type),
     15        parameter=md.control_type{i};
     16        %check this parameter is a field from model!
     17        if ~isfield(struct(md),parameter),
     18                error('displaysolutionparameters error message: one of the control type parameters does not exist!');
     19        end
     20        control_string=[control_string parameter ' and '];
     21end
     22control_string=control_string(1:length(control_string)-5);
     23disp(sprintf('         control_type: %s %s',control_string,'(list of parameters where inverse control is carried out; ex: {''drag''}, or {''drag'',''B''})'));
    1324disp(sprintf('         fit: (%i)      (''absolute: 0'', ''relative: 1'', or ''logarithmic: 2''. default is ''absolute: 0'', for each optimization steps)',length(md.fit)));
    1425disp(sprintf('         meanvel: %g      (velocity scaling factor when evaluating relative or logarithmic misfit)',md.meanvel));

  • issm/trunk/src/m/classes/public/displayqmu.m

    r1 r27  
    1010
    1111disp(sprintf('      ''%s''','qmu using Dakota'));
     12
    1213disp(sprintf('         variables:  (arrays of each variable class)'));
    1314fnames=fieldnames(md.variables);
    … …  
    1617        fnames{i},size(md.variables.(fnames{i})),class(md.variables.(fnames{i}))));
    1718end
     19
    1820disp(sprintf('         responses:  (arrays of each response class)'));
    1921fnames=fieldnames(md.responses);
    … …  
    2224        fnames{i},size(md.responses.(fnames{i})),class(md.responses.(fnames{i}))));
    2325end
     26
     27disp(sprintf('         qmu_method:  (array of dakota_method class)'));
     28for i=1:numel(md.qmu_method);
     29    if strcmp(class(md.qmu_method(i)),'dakota_method')
     30        disp(sprintf('            method%s :    ''%s''',...
     31            string_dim(md.qmu_method,i),md.qmu_method(i).method));
     32    end
     33end
     34
     35for i=1:numel(md.qmu_params)
     36    disp(sprintf('         qmu_params%s:  (array of method-independent parameters)',...
     37        string_dim(md.qmu_params,i)));
     38    fnames=fieldnames(md.qmu_params(i));
     39    maxlen=0;
     40    for j=1:numel(fnames)
     41        maxlen=max(maxlen,length(fnames{j}));
     42    end
     43   
     44    for j=1:numel(fnames)
     45        disp(sprintf(['            %-' num2str(maxlen+1) 's: %s'],...
     46            fnames{j},any2str(md.qmu_params(i).(fnames{j}))));
     47    end
     48end
     49
    2450disp(sprintf('         dakotaresults: 1x%i   {dvar,rfunc,scm,pcm,srcm,prcm}',length(md.dakotaresults)));
    2551if isempty(md.dakotain), disp(sprintf('         dakotain: N/A')); else disp(sprintf('         dakotain: not displayed (can be accessed by typing md.dakotain)'));end
    2652if isempty(md.dakotaout), disp(sprintf('         dakotaout: N/A')); else disp(sprintf('         dakotaout: not displayed (can be accessed by typing md.dakotaout)'));end
    2753if isempty(md.dakotadat), disp(sprintf('         dakotadat: N/A')); else disp(sprintf('         dakotadat: not displayed (can be accessed by typing md.dakotadat)'));end
    28 disp(sprintf('         method: ''%s'' (''local reliability'' or ''sampling'')',md.method));
    29 disp(sprintf('         seed: %g (random seed number)',md.seed));
    30 disp(sprintf('         samples: %i (number of samples)',md.samples));
    31 disp(sprintf('         sample_type: ''%s'' (type of sampling, ''random'' or ''lhs'')',md.sample_type));
    32 disp(sprintf('         interval_type: ''%s'' (FD ''forward'' or ''central'')',md.interval_type));
    33 disp(sprintf('         fd_gradient_step_size: %g (FD gradient step size, default= .001) ',md.fd_gradient_step_size));
    34 disp(sprintf('         evaluation_concurrency: %i (evaluation concurrency level)',md.evaluation_concurrency));
    35 disp(sprintf('         npart: %i (number of partitions for semi-descrete qmu)',md.npart));
    36 disp(sprintf('         analysis_driver    : ''%s'' (''matlab'' or name of driver for Dakota analysis)',md.analysis_driver));
    37 disp(sprintf('         analysis_components: ''%s'' (Matlab script file for Matlab driver)',md.analysis_components));
     54disp(sprintf('         npart   : %i (number of partitions for semi-descrete qmu)',md.npart));
     55disp(sprintf('         numrifts: %i (number of rifts for semi-descrete qmu)',md.numrifts));

  • issm/trunk/src/m/classes/public/extrude.m

    r1 r27  
    152152md.segmentonneumann_diag=segmentonneumann_diag;
    153153md.gridondirichlet_prog=project3d(md,md.gridondirichlet_prog,'node');
    154 md.dirichletvalues_prog=project3d(md,md.dirichletvalues_prog,'node',1);
     154md.dirichletvalues_prog=project3d(md,md.dirichletvalues_prog,'node');
    155155%md.segmentonneumann_prog=[tproj(md.segmentonneumann_prog(:,1)) tproj(md.segmentonneumann_prog(:,2)) tproj2d_el(md.segmentonneumann_prog(:,3))];
    156156%md.segmentonneumann_prog2=[tproj(md.segmentonneumann_prog2(:,1)) tproj(md.segmentonneumann_prog2(:,2)) tproj2d_el(md.segmentonneumann_prog2(:,3))];

  • issm/trunk/src/m/classes/public/ismodelselfconsistent.m

    r1 r27  
    55%      bool=ismodelselfconsistent(md,solutiontype,package),
    66
    7 if strcmpi(solutiontype,'mesh'),
    8         %this solution is a little special. It should come right after the md=model;  operation. So a lot less checks!
    9 
    10         bool=1;
    11         return;
    12 end
    13 
    147%tolerance we use in litmus tests for the consistency of the model
    158tolerance=10^-12;
    169if (nargin~=3  )
    17         IsModelSelfConsistentUsage();
    18         error(' ');
    19 end
    20 
     10        help ismodelselfconsistent
     11        error('ismodelselfconsistent error message: wrong usage');
     12end
     13
     14%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   COMMON CHECKS   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     15
     16%COUNTER
     17if md.counter<3,
     18        disp(['model ' md.name ' is not correctly configured. You forgot one step in the following sequence (mesh, geography, parameterize)!']);
     19        bool=0;return;
     20end
     21
     22%NAME
     23if isempty(md.name),
     24        disp(['model is not correctly configured: missing name!']);
     25        bool=0;return;
     26end
     27
     28%MESH
     29if md.numberofelements<=0,
     30        disp(['model ' md.name ' does not have any elements!']);
     31        bool=0; return;
     32end
     33if md.numberofgrids<=0,
     34        disp(['model ' md.name ' does not have any grids!']);
     35        bool=0; return;
     36end
     37
     38%ELEMENTSTYPE
     39if size(md.elements_type,1)~=md.numberofelements | size(md.elements_type,2)~=2,
     40        disp(['Types of elements have not been set properly, run setelementstype first'])
     41        bool=0;return;
     42end
     43if any(ones(md.numberofelements,1)-((md.elements_type(:,1)==hutterenum) + (md.elements_type(:,1)==macayealenum)  + (md.elements_type(:,1)==pattynenum)))
     44        disp(['Types of elements have not been set properly, run setelementstype first'])
     45        bool=0;return;
     46end
     47if any(ones(md.numberofelements,1)-((md.elements_type(:,2)==stokesenum) + (md.elements_type(:,2)==noneenum)))
     48        disp(['Types of elements have not been set properly, run setelementstype first'])
     49        bool=0;return;
     50end
     51if strcmpi(md.type,'2d'),
     52        if (ismember(pattynenum,md.elements_type(:,1)) |  ismember(stokesenum,md.elements_type(:,2))),
     53                disp(['For a 2d model, only MacAyeal''s and Hutter''s elements are allowed']);
     54                bool=0;return;
     55        end
     56end
     57
     58%NO NAN
     59fields={'numberofelements','numberofgrids','x','y','z','segmentonneumann_diag','drag','drag_type','p','q','gridondirichlet_diag',...
     60        'rho_ice','rho_water','B','elementoniceshelf','surface','thickness','bed','g','lowmem','sparsity','nsteps','maxiter',...
     61        'tolx','np','eps_rel','exclusive','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
     62for i=1:length(fields),
     63        if ~isempty(eval(['md.' char(fields(i))])),
     64                if find(isnan(eval(['md.' char(fields(i))]))),
     65                        disp(['model ' md.name ' has an NaN value in field ' char(fields(i)) '!']);
     66                        bool=0; return;
     67                end
     68        end
     69end
     70
     71%FIELDS > 0
     72fields={'numberofelements','numberofgrids','elements','segmentonneumann_diag','drag','drag_type','p','q','gridondirichlet_diag',...
     73        'rho_ice','rho_water','B','elementoniceshelf','thickness','g','eps_rel','eps_abs','nsteps','maxiter','tolx','exclusive',...
     74        'sparsity','lowmem','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
     75for i=1:length(fields),
     76        if ~isempty(eval(['md.' char(fields(i))])),
     77                if find((eval(['md.' char(fields(i))]))<0),
     78                        disp(['model ' md.name ' has a <0 value in field ' char(fields(i)) '!']);
     79                        bool=0; return;
     80                end
     81        end
     82end
     83if any(md.p<=0),
     84        disp(['model ' md.name ' has some p<0 friction coefficientsin sliding law']);
     85        bool=0; return;
     86end
     87
     88%FIELDS ~=0
     89fields={'numberofelements','numberofgrids','elements','segmentonneumann_diag','drag_type',...
     90        'rho_ice','rho_water','B','thickness','g','eps_rel','eps_abs','maxiter','tolx',...
     91        'sparsity','deltaH','DeltaH','timeacc','timedec'};
     92for i=1:length(fields),
     93        if ~isempty(eval(['md.' char(fields(i))])),
     94                if find((eval(['md.' char(fields(i))]))==0),
     95                        disp(['model ' md.name ' has a =0 value in field ' char(fields(i)) '!']);
     96                        bool=0; return;
     97                end
     98        end
     99end
     100
     101%SIZE NUMBEROFELEMENTS
     102fields={'elements','p','q','elementoniceshelf','n','elementonbed'};
     103for i=1:size(fields,2),
     104        if (size(eval(['md.' char(fields(i))]),1)~=md.numberofelements),
     105                disp(['model ' md.name ' field ' char(fields(i)) ' should be of size ' num2str(md.numberofelements) '!']);
     106                bool=0; return;
     107        end
     108end
     109
     110%SIZE NUMBEROFGRIDS
     111fields={'x','y','z','B','drag','gridondirichlet_diag','melting','accumulation','surface','thickness','bed','gridonbed','gridonsurface'};
     112for i=1:length(fields),
     113        if find(length(eval(['md.' char(fields(i))]))~=md.numberofgrids),
     114                disp(['model ' md.name ' field ' char(fields(i)) ' should be of size ' num2str(md.numberofgrids) '!']);
     115                bool=0; return;
     116        end
     117end
     118
     119%THICKNESS = SURFACE - BED
     120if find((md.thickness-md.surface+md.bed)>tolerance),
     121        disp(['model ' md.name ' violates the equality thickness=surface-bed!']);
     122        bool=0; return;
     123end
     124
     125%RIFTS
     126if md.numrifts,
     127        if ~strcmpi(md.type,'2d'),
     128                disp(['Models with rifts are only supported in 2d for now!']);
     129                bool=0;return;
     130        end
     131end
     132if ~isstruct(md.rifts),
     133        if ~isempty(find(md.segmentmarkers>=2)),
     134                %We have segments with rift markers, but no rift structure!
     135                disp(['model ' md.name ' should be processed for rifts (run meshprocessrifts)!']);
     136                bool=0; return;
     137        end
     138end
     139
     140%ARTIFICIAL DIFFUSIVITY
     141if ~isscalar(md.artificial_diffusivity),
     142        disp('artificial_diffusivity should be a scalar (1 or 0)');
     143        bool=0;return;
     144end
     145
     146%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  SOLUTION CHECKS  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     147
     148%DIAGNOSTIC
     149if strcmpi(solutiontype,'diagnostic')
     150
     151        %HUTTER ON ICESHELF WARNING
     152        if any(md.elements_type(:,1)==hutterenum & md.elementoniceshelf),
     153                disp(sprintf('\n !!! Warning: Hutter''s model is not consistent on ice shelves !!!\n'));
     154        end
     155
     156        %SINGULAR
     157        if ~any(md.gridondirichlet_diag),
     158                disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
     159                bool=0;return;
     160        end
     161
     162        %DIRICHLET VALUES
     163        if isempty(md.gridondirichlet_diag) | isempty(md.dirichletvalues_diag),
     164                disp(['model ' md.name ' is not well posed. Missing dirichlet values for diagnostic run']);
     165                bool=0;return;
     166        end
     167
     168        %DIRICHLET IF THICKNESS <= 0
     169        if ~isempty(find(md.thickness<=0)),
     170                pos=find(md.thickness<=0);
     171                if ~isempty(find(md.gridondirichlet_diag(pos)==0)),
     172                        disp(['model ' md.name ' has some grids with 0 thickness']);
     173                        bool=0; return;
     174                end
     175        end
     176end
     177
     178%PROGNOSTIC
     179if strcmp(solutiontype,'prognostic'),
     180        %old testing
     181        %       if isempty(find(md.gridondirichlet_diag)),
     182        %               disp(['model ' md.name ' diagnostic is not well posed (singular). You need at least one grid with fixed velocity!'])
     183        %               bool=0;return;
     184        %       end
     185        %       if isempty(find(md.gridondirichlet_prog)),
     186        %               disp(['model ' md.name ' prognostic is not well posed (singular). You need at least one grid with fixed thickness!'])
     187        %               bool=0;return;
     188        %       end
     189
     190        %VELOCITIES
     191        if (size(md.vx,1)~=md.numberofgrids | size(md.vy,1)~=md.numberofgrids),
     192                disp(['a 3d velocity is required. Run ''diagnostic'' solution first!'])
     193                bool=0; return;
     194        end
     195end
     196
     197%THERMAL TRANSIENT
     198if strcmp(solutiontype,'thermaltransient')
     199
     200        %INITIAL TEMPERATURE, MELTING AND ACCUMULATION
     201        if isempty(md.temperature),
     202                disp(['An  initial temperature is needed for a transient thermal computation'])
     203                bool=0;return;
     204        end
     205        if isstruct(md.temperature) |  isstruct(md.melting) | isstruct(md.accumulation),
     206                disp(['The initial temperature, melting or accumulation should be a list and not a structure'])
     207                bool=0;return;
     208        end
     209end
     210
     211%THERMAL STEADY AND THERMAL TRANSIENT
     212if strcmpi(solutiontype,'thermalsteady') |  strcmp(solutiontype,'thermaltransient'),
     213
     214        %EXTRUSION
     215        if strcmp(md.type,'2d'),
     216                disp(['For a ' solutiontype ' computation, the model must be 3d, extrude it first!'])
     217                bool=0;return;
     218        end
     219
     220        %VELOCITIES AND PRESSURE
     221        if (isempty(md.vx) | isnan(md.vx) | isempty(md.vy)  | isnan(md.vy) | isempty(md.vz) | isnan(md.vz)),
     222                disp(['a 3d velocity is required. Run ''diagnostic'' solution first!'])
     223                bool=0;return;
     224        end
     225        if (isempty(md.pressure) | isnan(md.pressure)),
     226                disp(['pressure is required. Run ''diagnostic'' solution first!'])
     227                bool=0;return;
     228        end
     229end
     230
     231%THERMAL TRANSIENT AND TRANSIENT
     232if strcmp(solutiontype,'thermaltransient') | strcmp(solutiontype,'transient'),
     233
     234        %DT and NDT
     235        fields={'dt','ndt'};
     236        for i=1:length(fields),
     237                if find((eval(['md.' char(fields(i))]))<0),
     238                        disp(['model ' md.name ' has a <0 value in field ' char(fields(i)) '!']);
     239                        bool=0; return;
     240                end
     241        end
     242
     243        %INITIAL TEMPERATURE
     244        if isstruct(md.temperature),
     245                disp(['The initial temperature should be empty or a list but not a structure'])
     246                bool=0;return;
     247        end
     248end
     249
     250%PARAMETERS
     251if strcmp(solutiontype,'parameters')
     252
     253        %OUTPUT
     254        if ~iscell(md.parameteroutput)
     255                disp(['parameteroutput field must be a cell, example {''strainrate'',''stress'',''deviatoricstress'',''viscousheating''}']);
     256                bool=0; return;
     257        end
     258        for i=1:length(md.parameteroutput)
     259                if ~strcmpi(md.parameteroutput(i),'strainrate') & ~strcmpi(md.parameteroutput(i),'stress')  & ~strcmpi(md.parameteroutput(i),'deviatoricstress') & ~strcmpi(md.parameteroutput(i),'viscousheating') ...
     260                                & ~strcmpi(md.parameteroutput(i),'pressure_elem') & ~strcmpi(md.parameteroutput(i),'stress_bed')  & ~strcmpi(md.parameteroutput(i),'stress_surface')
     261                        disp(['one of the parameteroutput is not supported yet']);
     262                        bool=0; return;
     263                end
     264        end
     265        %VELOCITY
     266        if ~(size(md.vx,1)==md.numberofgrids & size(md.vy,1)==md.numberofgrids & (size(md.vz,1)==md.numberofgrids | strcmpi(md.type,'2d')))
     267                disp(['velocities are required!']);
     268                bool=0;return;
     269        end
     270
     271        %ACCELERATION
     272        if strcmpi(md.type,'2d') & md.acceleration
     273                disp(['solution requires acceleration=0']);
     274                bool=0;return;
     275        end
     276
     277        %HUTTER
     278        if any(md.elements_type(:,1)==hutterenum);
     279                disp(['The model has Hutter''s elements. Impossible to compute parameters']);
     280                bool=0;return;
     281        end
     282end
     283
     284%CONTROL
     285if strcmpi(solutiontype,'control'),
     286
     287        %CONTROL TYPE
     288        if ~iscell(md.control_type),
     289                disp('control_type should be a cell array of strings');
     290                bool=0;return;
     291        end
     292
     293        %LENGTH CONTROL FIELDS
     294        if (length(md.maxiter)~=md.nsteps | length(md.optscal)~=md.nsteps | length(md.fit)~=md.nsteps)
     295                disp('maxiter, optscal and fit must have the length specified by nsteps')
     296                bool=0;return;
     297        end
     298
     299        %FIT
     300        if sum((double(md.fit==1) + double(md.fit==0) + double(md.fit==2))==1)~=md.nsteps
     301                disp('wrong fits: fit should be a vector of size nsteps holding 0, 1 and 2 only')
     302                bool=0;return;
     303        end
     304
     305        %OBSERVED VELOCITIES
     306        fields={'vx_obs','vy_obs'};
     307        for i=1:length(fields),
     308                if find(length(eval(['md.' char(fields(i))]))~=md.numberofgrids),
     309                        disp(['model ' md.name ' field ' char(fields(i)) ' should be of size ' num2str(md.numberofgrids) '!']);
     310                        bool=0; return;
     311                end
     312        end
     313end
     314
     315%QMU
    21316if strcmpi(solutiontype,'qmu'),
    22317        if ~strcmpi(md.cluster,'none'),
    … …  
    28323end
    29324
    30 if isempty(md.name),
    31         disp(['model is not correctly configured: missing name!']);
    32         bool=0;return;
    33 end
    34 
    35 if md.counter<3,
    36         disp(['model ' md.name ' is not correctly configured. You forgot one step in the following sequence (mesh, geography, parameterize)!']);
    37         bool=0;return;
    38 end
    39 
    40 if md.numberofelements<=0,
    41         disp(['model ' md.name ' does not have any elements!']);
    42         bool=0; return;
    43 end
    44 
    45 if md.numberofgrids<=0,
    46         disp(['model ' md.name ' does not have any grids!']);
    47         bool=0; return;
    48 end
    49 
    50 %check elementstype
    51 if size(md.elements_type,1)~=md.numberofelements | size(md.elements_type,2)~=2,
    52         disp(['Types of elements have not been set properly, run setelementstype first'])
    53         bool=0;return;
    54 end
    55 if any(ones(md.numberofelements,1)-((md.elements_type(:,1)==hutterenum) + (md.elements_type(:,1)==macayealenum)  + (md.elements_type(:,1)==pattynenum)))
    56         disp(['Types of elements have not been set properly, run setelementstype first'])
    57         bool=0;return;
    58 end
    59 if any(ones(md.numberofelements,1)-((md.elements_type(:,2)==stokesenum) + (md.elements_type(:,2)==noneenum)))
    60         disp(['Types of elements have not been set properly, run setelementstype first'])
    61         bool=0;return;
    62 end
    63 if strcmpi(md.type,'2d'),
    64         if (ismember(pattynenum,md.elements_type(:,1)) |  ismember(stokesenum,md.elements_type(:,2))),
    65                 disp(['For a 2d model, only MacAyeal''s and Hutter''s elements are allowed']);
    66                 bool=0;return;
    67         end
    68 end
    69 %check that if there are rifts, model is 2d
    70 if md.numrifts,
    71         if ~strcmpi(md.type,'2d'),
    72                 disp(['Models with rifts are only supported in 2d for now!']);
    73                 bool=0;return;
    74         end
    75 end
    76 
    77 %If there are rifts, check that the solver is set to mumps in parallel
    78 if md.numrifts,
    79         if ~strcmpi(md.cluster,'none')
    80                 if isempty(findstr('aijmumps',md.solverstring)),
    81                         disp(['For a 2d model with rifts, running on a cluster, you should use a direct solver like MUMPS!']);
     325%MESH
     326if strcmpi(solutiontype,'mesh'),
     327        %this solution is a little special. It should come right after the md=model;  operation. So a lot less checks!
     328
     329        bool=1;
     330        return;
     331end
     332
     333%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  PACKAGE CHECKS   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     334
     335%CIELO
     336if strcmpi(package,'cielo'),
     337
     338        %NAN VALUES
     339        fields={'time','sparsity'};
     340        for i=1:length(fields),
     341                if ~isempty(eval(['md.' char(fields(i))])),
     342                        if find(isnan(eval(['md.' char(fields(i))]))),
     343                                disp(['model ' md.name ' has an NaN value in field ' char(fields(i)) '!']);
     344                                bool=0; return;
     345                        end
     346                end
     347        end
     348
     349        %FIELD > 0
     350        fields={'time','sparsity'};
     351        for i=1:length(fields),
     352                if ~isempty(eval(['md.' char(fields(i))])),
     353                        if find((eval(['md.' char(fields(i))]))<0),
     354                                disp(['model ' md.name ' has a <0 value in field ' char(fields(i)) '!']);
     355                                bool=0; return;
     356                        end
     357                end
     358        end
     359
     360        %FIELD ~= 0
     361        fields={'time','sparsity'};
     362        for i=1:length(fields),
     363                if ~isempty(eval(['md.' char(fields(i))])),
     364                        if find((eval(['md.' char(fields(i))]))==0),
     365                                disp(['model ' md.name ' has a =0 value in field ' char(fields(i)) '!']);
     366                                bool=0; return;
     367                        end
     368                end
     369        end
     370
     371        %SPARSITY BETWEEN 0 AND 1
     372        if ( (md.sparsity<=0) | (md.sparsity>1)),
     373                disp(['model ' md.name ' sparsity should be inside the ]0 1] range']);
     374                bool=0; return;
     375        end
     376
     377        %RIFTS
     378        if md.numrifts,
     379                if ~strcmpi(md.cluster,'none')
     380                        if isempty(findstr('aijmumps',md.solverstring)),
     381                                disp(['For a 2d model with rifts, running on a cluster, you should use a direct solver like MUMPS!']);
     382                                bool=0;return;
     383                        end
     384                end
     385        end
     386
     387        %CONNECTIVITY
     388        if strcmpi(md.type,'2d'),
     389                if md.connectivity<9,
     390                        disp('connectivity should be at least 9 for 2d models');
    82391                        bool=0;return;
    83392                end
    84393        end
    85 end
    86 
    87 %check elementstype in diagnostic
    88 if strcmpi(solutiontype,'diagnostic')
    89         if any(md.elements_type(:,1)==hutterenum & md.elementoniceshelf),
    90                 disp(sprintf('\n !!! Warning: Hutter''s model is not consistent on ice shelves !!!\n'));
    91         end
    92 end
    93 
    94 %Check that mesh is 3d
    95 if strcmpi(solutiontype,'thermalsteady') |  strcmp(solutiontype,'thermaltransient'),
    96         if strcmp(md.type,'2d'),
    97                 disp(['For a ' solutiontype ' computation, the model must be 3d, extrude it first!'])
    98                 bool=0;return;
    99         end
    100 end
    101 
    102 %check artificial_diffusivity is an integer
    103 if ~isscalar(md.artificial_diffusivity),
    104         disp('artificial_diffusivity should be a scalar (1 or 0)');
    105         bool=0;return;
    106 end
    107 
    108 %check on connectivity
    109 if strcmpi(md.type,'2d'),
    110         if md.connectivity<9,
    111                 disp('connectivity should be at least 9 for 2d models');
    112                 bool=0;return;
    113         end
    114 end
    115 
    116 if strcmpi(md.type,'3d'),
    117         if md.connectivity<24,
    118                 disp('connectivity should be at least 24 for 3d models');
    119                 bool=0;return;
    120         end
    121 end
    122 
    123 
    124 
    125 
    126 
    127 %Check previous computation
    128 if strcmp(solutiontype,'thermaltransient') | strcmp(solutiontype,'thermalsteady')
    129         if (isempty(md.vx) | isempty(md.vy) | isempty(md.vz)),
    130                 disp(['a 3d velocity is required. Run ''diagnostic'' solution first!'])
    131                 bool=0;return;
    132         end
    133 end
    134 
    135 if strcmp(solutiontype,'thermaltransient') | strcmp(solutiontype,'thermalsteady')
    136         if (isempty(md.pressure) | isnan(md.pressure)),
    137                 disp(['pressure is required. Run ''diagnostic'' solution first!'])
    138                 bool=0;return;
    139         end
    140        
    141         if (isempty(md.vx) | isnan(md.vx) | isempty(md.vy)  | isnan(md.vy) | isempty(md.vz) | isnan(md.vz)),
    142                 disp(['velocities are required!']);
    143                 bool=0;return;
    144         end
    145 end
    146 
    147 if strcmp(solutiontype,'thermaltransient')
    148         if isempty(md.temperature),
    149                 disp(['An  initial temperature is needed for a transient thermal computation'])
    150                 bool=0;return;
    151         end
    152         if isstruct(md.temperature) |  isstruct(md.melting),
    153                 disp(['The initial temperature or melting should be a list and not a structure'])
    154                 bool=0;return;
    155         end
    156 
    157 end
    158 
    159 if strcmp(solutiontype,'transient')
    160         if isstruct(md.temperature),
    161                 disp(['The initial temperature should be empty or a list but not a structure'])
    162                 bool=0;return;
    163         end
    164 
    165 end
    166 
    167 if strcmp(solutiontype,'parameters')
    168         if ~(size(md.vx,1)==md.numberofgrids & size(md.vy,1)==md.numberofgrids & (size(md.vz,1)==md.numberofgrids | strcmpi(md.type,'2d')))
    169                 disp(['velocities are required!']);
    170                 bool=0;return;
    171         end
    172         if strcmpi(md.type,'2d') & md.acceleration
    173                 disp(['solution requires acceleration=0']);
    174                 bool=0;return;
    175         end
    176         if any(md.elements_type(:,1)==hutterenum);
    177                 disp(['The model has Hutter''s elements. Impossible to compute parameters']);
    178                 bool=0;return;
    179         end
    180 end
    181 
    182 
    183 %Check that problem is not singular
    184 if strcmp(solutiontype,'diagnostic'),
    185         if isempty(find(md.gridondirichlet_diag)),
    186                 disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
    187                 bool=0;return;
    188         end
    189 end
    190 
    191 %Check we have values for dirichletvalues_diag
    192 if strcmp(solutiontype,'diagnostic'),
    193         if isempty(md.gridondirichlet_diag) | isempty(md.dirichletvalues_diag),
    194                 disp(['model ' md.name ' is not well posed. Missing dirichlet values for diagnostic run']);
    195                 bool=0;return;
    196         end
    197 end
    198 
    199 
    200 
    201 if strcmp(solutiontype,'prognostic'),
    202         %old testing
    203 %       if isempty(find(md.gridondirichlet_diag)),
    204 %               disp(['model ' md.name ' diagnostic is not well posed (singular). You need at least one grid with fixed velocity!'])
    205 %               bool=0;return;
    206 %       end
    207 %       if isempty(find(md.gridondirichlet_prog)),
    208 %               disp(['model ' md.name ' prognostic is not well posed (singular). You need at least one grid with fixed thickness!'])
    209 %               bool=0;return;
    210 %       end
    211         if (size(md.vx,1)~=md.numberofgrids | size(md.vy,1)~=md.numberofgrids),
    212                 disp(['a 3d velocity is required. Run ''diagnostic'' solution first!'])
    213                 bool=0; return;
    214         end
    215 end
    216 %some checks for control methods
    217 if strcmp(solutiontype,'control')
    218         %check fields used by control
    219         if (length(md.maxiter)~=md.nsteps | length(md.optscal)~=md.nsteps | length(md.fit)~=md.nsteps)
    220                 disp('maxiter, optscal and fit must have the length specified by nsteps')
    221                 bool=0;return;
    222         end
    223         if sum((double(md.fit==1) + double(md.fit==0) + double(md.fit==2))==1)~=md.nsteps
    224                 disp('wrong fits: fit should be a vector of size nsteps holding 0, 1 and 2 only')
    225                 bool=0;return;
    226         end
    227 end
    228 
    229 
    230 %Check that no NaN values popped up:
    231 fields={'numberofelements','numberofgrids','x','y','z','segmentonneumann_diag','drag','drag_type','p','q','gridondirichlet_diag',...
    232         'rho_ice','rho_water','B','elementoniceshelf','surface','thickness','bed','g','lowmem','sparsity','nsteps','maxiter',...
    233                 'tolx','np','time','eps_rel','exclusive','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
    234 for i=1:length(fields),
    235         if ~isempty(eval(['md.' char(fields(i))])),
    236                 if find(isnan(eval(['md.' char(fields(i))]))),
    237                         disp(['model ' md.name ' has an NaN value in field ' char(fields(i)) '!']);
    238                         bool=0; return;
    239                 end
    240         end
    241 end
    242 
    243 %Check that these fields are > 0
    244 fields={'numberofelements','numberofgrids','elements','segmentonneumann_diag','drag','drag_type','p','q','gridondirichlet_diag',...
    245          'rho_ice','rho_water','B','elementoniceshelf','thickness','g','eps_rel','eps_abs','nsteps','maxiter','tolx','np','time','exclusive',...
    246                  'sparsity','lowmem','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
    247 for i=1:length(fields),
    248         if ~isempty(eval(['md.' char(fields(i))])),
    249                 if find((eval(['md.' char(fields(i))]))<0),
    250                         disp(['model ' md.name ' has a <0 value in field ' char(fields(i)) '!']);
    251                         bool=0; return;
    252                 end
    253         end
    254 end
    255 
    256 %Check that these fields are ~=0
    257 fields={'numberofelements','numberofgrids','elements','segmentonneumann_diag','drag_type',...
    258          'rho_ice','rho_water','B','thickness','g','eps_rel','eps_abs','maxiter','tolx','np','time',...
    259                  'sparsity','deltaH','DeltaH','timeacc','timedec'};
    260 for i=1:length(fields),
    261         if ~isempty(eval(['md.' char(fields(i))])),
    262                 if find((eval(['md.' char(fields(i))]))==0),
    263                         disp(['model ' md.name ' has a =0 value in field ' char(fields(i)) '!']);
    264                         bool=0; return;
    265                 end
    266         end
    267 end
    268 
    269 %Check that these fields are of size md.numberofelements
    270 fields={'elements','p','q','elementoniceshelf','n','elementonbed'};
    271 for i=1:size(fields,2),
    272         if (size(eval(['md.' char(fields(i))]),1)~=md.numberofelements),
    273                 disp(['model ' md.name ' field ' char(fields(i)) ' should be of size ' num2str(md.numberofelements) '!']);
    274                 bool=0; return;
    275         end
    276 end
    277 
    278 %Check that these fields are of size md.numberofgrids
    279 fields={'melting','x','y','z','B','drag','gridondirichlet_diag','vx_obs','vy_obs','surface','thickness','bed','gridonbed','gridonsurface'};
    280 for i=1:length(fields),
    281         if find(length(eval(['md.' char(fields(i))]))~=md.numberofgrids),
    282                 disp(['model ' md.name ' field ' char(fields(i)) ' should be of size ' num2str(md.numberofgrids) '!']);
    283                 bool=0; return;
    284         end
    285 end
    286 
    287 %Make sure lowmem is either 0 or 1
    288 if ((md.lowmem ~= 1) & (md.lowmem~=0)),
    289         disp(['model ' md.name ' lowmem field should be 0 or 1']);
    290         bool=0; return;
    291 end
    292 
    293 %Make sure sparsity is between 0 and 1:
    294 if ( (md.sparsity<=0) | (md.sparsity>1)),
    295         disp(['model ' md.name ' sparsity should be inside the ]0 1] range']);
    296         bool=0; return;
    297 end
    298 
    299 %Make sure thickness=surface-bed;
    300 if find((md.thickness-md.surface+md.bed)>tolerance),
    301         disp(['model ' md.name ' violates the equality thickness=surface-bed!']);
    302         bool=0; return;
    303 end
    304 
    305 %Check that rifts were processed correctly:
    306 if ~isstruct(md.rifts),
    307         if ~isempty(find(md.segmentmarkers>=2)),
    308                 %We have segments with rift markers, but no rift structure!
    309                 disp(['model ' md.name ' should be processed for rifts (run meshprocessrifts)!']);
    310                 bool=0; return;
    311         end
    312 end
    313 
    314 
    315 %Check that thickness is not null. If is it, then the corresponding grids should be spc'd.
    316 if ~isempty(find(md.thickness<=0)),
    317         pos=find(md.thickness<=0);
    318         if ~isempty(find(md.gridondirichlet_diag(pos)==0)),
    319                 disp(['model ' md.name ' has some grids with 0 thickness']);
    320                 bool=0; return;
    321         end
    322 end
    323 
    324 %Check that p>0 in the sliding law u = k * sigma ^ p * Neff ^ q :
    325 if ~isempty(find(md.p<=0)),
    326         disp(['model ' md.name ' has some p<0 friction coefficientsin sliding law']);
    327         bool=0; return;
    328 end
    329 
    330 
    331 %Check parameteroutput
    332 if ~iscell(md.parameteroutput)
    333         disp(['parameteroutput field must be a cell, example {''strainrate'',''stress'',''deviatoricstress'',''viscousheating''}']);
    334         bool=0; return;
    335 end
    336 
    337 
    338 for i=1:length(md.parameteroutput)
    339         if ~strcmpi(md.parameteroutput(i),'strainrate') & ~strcmpi(md.parameteroutput(i),'stress')  & ~strcmpi(md.parameteroutput(i),'deviatoricstress') & ~strcmpi(md.parameteroutput(i),'viscousheating') ...
    340                  & ~strcmpi(md.parameteroutput(i),'pressure_elem') & ~strcmpi(md.parameteroutput(i),'stress_bed')  & ~strcmpi(md.parameteroutput(i),'stress_surface')
    341                 disp(['one of the parameteroutput is not supported yet']);
    342                 bool=0; return;
    343         end
    344 end
    345 
    346 %check parallel settings
    347 if ~strcmpi(md.cluster,'none'),
    348         if ~strcmpi(package,'cielo'),
    349                 disp(['parallel runs only allowed for ''cielo'' package']);
    350                 bool=0; return;
    351         end
    352 end
    353 
    354 %Put as many checks as you want here:
    355 
     394        if strcmpi(md.type,'3d'),
     395                if md.connectivity<24,
     396                        disp('connectivity should be at least 24 for 3d models');
     397                        bool=0;return;
     398                end
     399        end
     400
     401        %NP
     402        if md.np==0,
     403                disp(['model ' md.name ' has a =0 value in field ' md.np ' !']);
     404        elseif md.np<0,
     405                disp(['model ' md.name ' has a negative value in field ' md.np ' !']);
     406        end
     407
     408        %LOWMEM = 0 or 1
     409        if ((md.lowmem ~= 1) & (md.lowmem~=0)),
     410                disp(['model ' md.name ' lowmem field should be 0 or 1']);
     411                bool=0; return;
     412        end
     413end
    356414
    357415%No problems, just return 1;
    358416bool=1;
    359417return;
    360 
    361 function IsModelSelfConsistentUsage()
    362 disp(' ');
    363 disp('   IsModelSelfConsistent usage: flag=IsModelSelfConsistent(model)');
    364 disp('         where model is an instance of the @model class, and flag a boolean');
    365 
    366 
    367 
    368 

  • issm/trunk/src/m/classes/public/isresultconsistent.m

    r1 r27  
    2121                if ~isnan(md.penalties),
    2222                        for i=1:size(md.penalties,1)
    23                                 for n=2:md.numlayers
     23                                for n=2:size(md.penalties,2)
    2424                                        relativex=(md.vx(md.penalties(i,1))-md.vx(md.penalties(i,n)))./(md.vx(md.penalties(i,1)));
    2525                                        if md.vx(md.penalties(i,1))==md.vx(md.penalties(i,n)), relativex=0; end
    … …  
    7474                if ~isnan(md.penalties),
    7575                        for i=1:size(md.penalties,1)
    76                                 for n=2:md.numlayers
     76                                for n=2:size(md.penalties,2)
    7777                                        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)));
    7878                                        if md.transient_results(iter).vx(md.penalties(i,1))==md.transient_results(iter).vx(md.penalties(i,n)), relativex=0; end

  • issm/trunk/src/m/classes/public/parameterize.m

    r1 r27  
    2727end
    2828
    29 %Try and run parameter file. We try to capture the error message, but this is not allowed for
    30 %lower versions of matlab.
     29%Try and run parameter file.
     30eval(['system( '' cp ' parametername  ' TemporaryParameterFile.m'');']);
    3131
    32 stringfile=readfile(parametername);
    33 stringlines=strsplit(stringfile,char(10)); %split into lines
    34 
    35 for i=1:length(stringlines),
    36         a=stringlines{i};
    37         try,
    38                 eval(stringlines{i});
    39         catch me,
    40                 disp(' ');
    41                 disp(['Error in ==> ' parametername ' at ' num2str(i)]);
    42                 disp(me.message);
    43                 disp(stringlines{i});
    44                 error(' ');
     32try,
     33        TemporaryParameterFile
     34        system('rm TemporaryParameterFile.m');
     35catch me,
     36        system('rm TemporaryParameterFile.m');
     37        me2=struct('message',me.message,'stack',me.stack);
     38        for i=1:length(me2.stack)
     39                if (length(me2.stack(i).file)>=24 & strcmp(me2.stack(i).file(end-23:end),'TemporaryParameterFile.m'))
     40                        me2.stack(i).file=[me2.stack(i).file(1:end-24) parametername];
     41                end
     42                if strcmp(me2.stack(i).name,'TemporaryParameterFile'),
     43                        me2.stack(i).name=parametername;
     44                end
    4545        end
     46        rethrow(me2);
    4647end
    4748

  • issm/trunk/src/m/classes/public/plot/applyoptions.m

    r1 r27  
    4444if ~isnan(options_structure.view),
    4545        view(options_structure.view);
     46else
     47        if strcmpi(md.type,'3d') & isnan(options_structure.layer),
     48                view(3);
     49        else
     50                view(2);
     51        end
    4652end
    4753
    … …  
    162168end
    163169
     170%streamliness
     171if iscell(options_structure.streamlines) | ~isnan(options_structure.streamlines),
     172        plot_streamlines(md,options_structure);
     173end
     174
    164175%contours
    165176if iscell(options_structure.contourlevels) | ~isnan(options_structure.contourlevels),

  • issm/trunk/src/m/classes/public/plot/imagescnan.m

    r1 r27  
    8080
    8181%   Copyright 2008 Carlos Adrian Vargas Aguilera
    82 %   $Revision: 1.1 $  $Date: 2009/03/24 21:05:20 $
     82%   $Revision: 1.1 $  $Date: 2009/04/03 22:56:05 $
    8383
    8484%   Written by

  • issm/trunk/src/m/classes/public/plot/parse_options.m

    r1 r27  
    1717end
    1818
     19%density
     20densityvalues=findarg(optionstring,'density');
     21if ~isempty(densityvalues),
     22        options_struct.density=abs(ceil(densityvalues(1).value));
     23else
     24        options_struct.density=NaN;
     25end
     26
    1927%Section profile
    2028sectionvalues=findarg(optionstring,'sectionvalue');
    … …  
    3846
    3947%Show section
    40 showsections=findarg(optionstring,'showsection');
    41 if ~isempty(showsections),
    42         if ischar(showsections(1).value),
    43                 options_struct.showsection=showsections(1).value;
    44         else
    45                 options_struct.showsection=0;
     48showsection=findarg(optionstring,'showsection');
     49if ~isempty(showsection),
     50        if ischar(showsection(1).value) & strcmpi(showsection(1).value,'on'),
     51                options_struct.showsection=4;
     52        else
     53                options_struct.showsection=showsection(1).value;
    4654        end
    4755else
    … …  
    115123else
    116124        options_struct.smooth=NaN;
     125end
     126
     127%streamlines
     128streamlines_values=findarg(optionstring,'streamlines');
     129if ~isempty(streamlines_values),
     130        options_struct.streamlines=streamlines_values.value;
     131else
     132        options_struct.streamlines=NaN;
    117133end
    118134
    … …  
    226242        options_struct.view=viewvalues.value;
    227243else
    228         if strcmpi(md.type,'3d') & isnan(options_struct.layer),
    229                 options_struct.view=3;
    230         else
    231                 options_struct.view=2;
    232         end
     244        options_struct.view=NaN;
    233245end
    234246

  • issm/trunk/src/m/classes/public/plot/plot_basaldrag.m

    r1 r27  
    1 function plot_basaldrag(md,options_structure,width,i);
     1function plot_basaldrag(md,options_structure,width,i,type);
    22%PLOT_BASALDRAG - plot basal drag
    33%
    44%   Usage:
    5 %      plot_basaldrag(md,options_structure,width,i);
     5%      plot_basaldrag(md,options_structure,width,i,type);
    66%
    77%   See also: PLOTMODEL
    … …  
    2020
    2121%compute horizontal velocity
    22 ub=sqrt(md.vx.^2+md.vy.^2)/md.yts;
     22if strcmpi(type,'basal_drag')
     23        ub=sqrt(md.vx.^2+md.vy.^2)/md.yts;
     24elseif strcmpi(type,'basal_dragx')
     25        ub=md.vx/md.yts;
     26elseif strcmpi(type,'basal_dragy')
     27        ub=md.vy/md.yts;
     28end
    2329
    2430%compute basal drag
    2531drag=(md.g*(md.rho_ice*md.thickness+md.rho_water*md.bed)).^r.*(md.drag).^2.*ub.^s/1000;
    2632
    27 %process data and model
    28 [x y z elements is2d]=processmesh(md,options_structure);
    29 [basal_drag isongrid]=processdata(md,drag,options_structure);
     33%Figure out if this is a Section plot
     34if ~isnan(options_structure.sectionvalue)
     35        plot_section(md,drag,options_structure,width,i);
     36        return;
     37else
    3038
    31 %plot mesh quivervel
    32 subplot(width,width,i);
     39        %process data and model
     40        [x y z elements is2d]=processmesh(md,options_structure);
     41        [basal_drag isongrid]=processdata(md,drag,options_structure);
    3342
    34 %edgecolor?
    35 if ~isnan(options_structure.edgecolor),
    36         edgecolor=options_structure.edgecolor;
    37 else
    38         edgecolor='none';
     43        %plot basaldrag
     44        subplot(width,width,i);
     45        plot_unit(x,y,z,elements,basal_drag,isongrid,is2d,options_structure);
     46
     47        %apply options
     48        if isnan(options_structure.title)
     49                options_structure.title='Basal drag [kPa]';
     50        end
     51        if isnan(options_structure.view)
     52                options_structure.view=2;
     53        end
     54        applyoptions(md,basal_drag,options_structure);
     55
    3956end
    40 
    41 %plot basal basal drag
    42 A=elements(:,1); B=elements(:,2); C=elements(:,3);
    43 
    44 patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData',basal_drag,'FaceColor','interp','EdgeColor',edgecolor);
    45 
    46 %apply options
    47 if isnan(options_structure.title)
    48         options_structure.title='Basal drag [kPa]';
    49 end
    50 options_structure.view=2;
    51 applyoptions(md,basal_drag,options_structure);

  • issm/trunk/src/m/classes/public/plot/plot_boundaries.m

    r1 r27  
    4141        options_structure.colorbar=0;
    4242end
     43if isnan(options_structure.view)
     44        options_structure.view=2;
     45end
    4346applyoptions(md,[],options_structure);

  • issm/trunk/src/m/classes/public/plot/plot_contour.m

    r1 r27  
    4646
    4747%initialization of some variables
    48 numberofelements=md.numberofelements;
     48numberofelements=size(index,1);
    4949elementslist=1:numberofelements;
    5050c=[];
    … …  
    231231
    232232%labels?
    233 if ~strcmpi(options_structure.contourticks,'off')
     233if (~strcmpi(options_structure.contourticks,'off') & ~isempty(c) & ~isempty(h))
    234234        if ~isnan(options_structure.contourcolor)
    235235                clabel(c,h,'color',color);

  • issm/trunk/src/m/classes/public/plot/plot_drivingstress.m

    r1 r27  
    55%      plot_drivingstress(md,options_structure,width,i);
    66%
    7 %   See also: PLOTMODEL
     7%   See also: PLOTMODEL, PLOT_UNIT, PLOT_DEALER
    88
    99%get driving stress
    … …  
    1717%plot mesh quivervel
    1818subplot(width,width,i);
    19 
    20 %edgecolor?
    21 if ~isnan(options_structure.edgecolor),
    22         edgecolor=options_structure.edgecolor;
    23 else
    24         edgecolor='none';
    25 end
    26 
    27 %element data
    28 if ~isongrid
    29         if is2d
    30                 A=elements(:,1); B=elements(:,2); C=elements(:,3);
    31                 patch( 'Faces', [A B C], 'Vertices', [x y z],'CData', dstress,'FaceColor','flat','EdgeColor',edgecolor);
    32         else
    33                 A=elements(:,1); B=elements(:,2); C=elements(:,3); D=elements(:,4); E=elements(:,5); F=elements(:,6);
    34                 patch( 'Faces', [A B C], 'Vertices', [x y z],'CData', dstress,'FaceColor','flat','EdgeColor',edgecolor);
    35                 patch( 'Faces', [D E F], 'Vertices', [x y z],'CData', dstress,'FaceColor','flat','EdgeColor',edgecolor);
    36                 patch( 'Faces', [A B E D], 'Vertices', [x y z],'CData', dstress,'FaceColor','flat','EdgeColor',edgecolor);
    37                 patch( 'Faces', [B E F C ], 'Vertices', [x y z],'CData', dstress,'FaceColor','flat','EdgeColor',edgecolor);
    38                 patch( 'Faces', [C A D F ], 'Vertices', [x y z],'CData', dstress,'FaceColor','flat','EdgeColor',edgecolor);
    39         end
    40 %grid data
    41 elseif isongrid
    42         if is2d
    43                 A=elements(:,1); B=elements(:,2); C=elements(:,3);
    44                 patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', dstress,'FaceColor','interp','EdgeColor',edgecolor);
    45         else
    46                 if options_structure.layer>=1,
    47                         A=elements(:,1); B=elements(:,2); C=elements(:,3);
    48                         patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', dstress,'FaceColor','interp','EdgeColor',edgecolor);
    49                 else
    50                         A=elements(:,1); B=elements(:,2); C=elements(:,3); D=elements(:,4); E=elements(:,5); F=elements(:,6);
    51                         patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', dstress,'FaceColor','interp','EdgeColor',edgecolor);
    52                         patch( 'Faces', [D E F], 'Vertices', [x y z],'FaceVertexCData', dstress,'FaceColor','interp','EdgeColor',edgecolor);
    53                         patch( 'Faces', [A B E D], 'Vertices', [x y z],'FaceVertexCData', dstress,'FaceColor','interp','EdgeColor',edgecolor);
    54                         patch( 'Faces', [B E F C ], 'Vertices', [x y z],'FaceVertexCData', dstress,'FaceColor','interp','EdgeColor',edgecolor);
    55                         patch( 'Faces', [C A D F ], 'Vertices', [x y z],'FaceVertexCData', dstress,'FaceColor','interp','EdgeColor',edgecolor);
    56                 end
    57         end
    58 end
     19plot_unit(x,y,z,elements,dstress,isongrid,is2d,options_structure)
    5920
    6021%apply options
    … …  
    6223        options_structure.title='Driving stress [kPa]';
    6324end
     25if isnan(options_structure.view)
     26        options_structure.view=2;
     27end
    6428applyoptions(md,dstress,options_structure);

  • issm/trunk/src/m/classes/public/plot/plot_elementnumbering.m

    r1 r27  
    55%      plot_elementnumbering(md,options_structure,width,i);
    66%
    7 %   See also: PLOTMODEL
     7%   See also: PLOTMODEL, PLOT_UNIT, PLOT_DEALER
    88
    99subplot(width,width,i);

  • issm/trunk/src/m/classes/public/plot/plot_elementstype.m

    r1 r27  
    7676        if ~isempty(posS)
    7777                A=elements(posS,1); B=elements(posS,2); C=elements(posS,3); D=elements(posS,4); E=elements(posS,5); F=elements(posS,6);
    78         %       p4=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData', stokesenum,'FaceColor','flat','EdgeColor',edgecolor,'EdgeAlpha',alpha,'FaceAlpha',alpha);
    79         %       patch( 'Faces', [D E F], 'Vertices', [x y z],'CData', stokesenum,'FaceColor','flat','EdgeColor',edgecolor,'EdgeAlpha',alpha,'FaceAlpha',alpha);
    80         %       patch( 'Faces', [A B E D], 'Vertices', [x y z],'CData', stokesenum,'FaceColor','flat','EdgeColor',edgecolor,'EdgeAlpha',alpha,'FaceAlpha',alpha);
    81         %       patch( 'Faces', [B E F C ], 'Vertices', [x y z],'CData', stokesenum,'FaceColor','flat','EdgeColor',edgecolor,'EdgeAlpha',alpha,'FaceAlpha',alpha);
    82         %       patch( 'Faces', [C A D F ], 'Vertices', [x y z],'CData', stokesenum,'FaceColor','flat','EdgeColor',edgecolor,'EdgeAlpha',alpha,'FaceAlpha',alpha);
    83         %       legend([p1 p2 p3 p4],'Hutter''s elements','MacAyeal''s elements','Pattyn''s elements','Stokes''s elements');
     78                p4=patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceColor','none','EdgeColor','black');
     79                patch( 'Faces', [D E F], 'Vertices', [x y z],'FaceColor','none','EdgeColor','black');
     80                patch( 'Faces', [A B E D], 'Vertices', [x y z],'FaceColor','none','EdgeColor','black');
     81                patch( 'Faces', [B E F C ], 'Vertices', [x y z],'FaceColor','none','EdgeColor','black');
     82                patch( 'Faces', [C A D F ], 'Vertices', [x y z],'FaceColor','none','EdgeColor','black');
     83                legend([p1 p2 p3 p4],'Hutter''s elements','MacAyeal''s elements','Pattyn''s elements','Stokes''s elements');
    8484        else
    8585                legend([p1 p2 p3],'Hutter''s elements','MacAyeal''s elements','Pattyn''s elements');

  • issm/trunk/src/m/classes/public/plot/plot_penalties.m

    r1 r27  
    2222if ~strcmpi(md.type,'3d'),
    2323        error('no penalties to plot for ''2d'' model');
    24 elseif isnan(md.penalties) | isempty(md.penalties)
     24elseif isempty(md.penalties),
    2525        disp('no penalty applied in this model');
    2626        return;

  • issm/trunk/src/m/classes/public/plot/plot_quiver.m

    r1 r27  
    1616                                       
    1717%quiver 2d
     18if ~isnan(options_structure.density)
     19        x=x(1:options_structure.density:end);
     20        y=y(1:options_structure.density:end);
     21        vx=vx(1:options_structure.density:end);
     22        vy=vy(1:options_structure.density:end);
     23end
    1824quiver(x,y,vx,vy);
    1925

  • issm/trunk/src/m/classes/public/plot/plot_quiver3.m

    r1 r27  
    1717
    1818%quiver 3d
     19if ~isnan(options_structure.density)
     20        x=x(1:options_structure.density:end);
     21        y=y(1:options_structure.density:end);
     22        z=z(1:options_structure.density:end);
     23        vx=vx(1:options_structure.density:end);
     24        vy=vy(1:options_structure.density:end);
     25        vz=vz(1:options_structure.density:end);
     26end
    1927quiver3(x,y,z,vx,vy,vz);
    2028

  • issm/trunk/src/m/classes/public/plot/plot_quivervel.m

    r1 r27  
    55%      plot_quivervel(md,options_structure,width,i);
    66%
    7 %   See also: PLOTMODEL
     7%   See also: PLOTMODEL, PLOT_UNIT, PLOT_DEALER
    88
    99%process data and model
    … …  
    1515end
    1616
    17 %edgecolor?
    18 if ~isnan(options_structure.edgecolor),
    19         edgecolor=options_structure.edgecolor;
    20 else
    21         edgecolor='none';
     17%density
     18if ~isnan(options_structure.density)
     19        x=x(1:options_structure.density:end);
     20        y=y(1:options_structure.density:end);
     21        z=z(1:options_structure.density:end);
     22        vx=vx(1:options_structure.density:end);
     23        vy=vy(1:options_structure.density:end);
     24        vz=vz(1:options_structure.density:end);
    2225end
    2326
    2427%plot
    2528subplot(width,width,i);
     29colormap('default')
     30plot_unit(x,y,z,elements,sqrt(vx.^2+vy.^2),isongrid,is2d,options_structure)
    2631
    2732if is2d
    2833        hold on
    29         colormap('default')
    30         A=elements(:,1); B=elements(:,2); C=elements(:,3);
    31         patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor);
    3234        colorbar;
    3335        quiver(x,y,vx,vy,'k')
    … …  
    3537else
    3638        hold on
    37         colormap('default')
    38         A=elements(:,1); B=elements(:,2); C=elements(:,3); D=elements(:,4); E=elements(:,5); F=elements(:,6);
    39         patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor);
    40         patch( 'Faces', [D E F], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor);
    41         patch( 'Faces', [A B E D], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor);
    42         patch( 'Faces', [B E F C ], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor);
    43         patch( 'Faces', [C A D F ], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor);
    4439        colorbar;
    4540        quiver3(x,y,z,vx,vy,vz);

  • issm/trunk/src/m/classes/public/plot/plot_section.m

    r1 r27  
    88
    99%How many subplots?
    10 if strcmpi(options_structure.showsection,'yes')
     10if ~isnan(options_structure.showsection)
    1111
    1212        %Compute the indexes of the 2 plots (one for the sectionvalue and one for showsection
    … …  
    2020end
    2121
    22 %check on arguments
    23 if (iscell(data) | isempty(data)),
    24         error('plot error message: data provided is empty');
    25 end
    26 if md.numberofgrids==size(md.elements,1),
    27         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')
    28 end
    2922
    30 %smoothing?
    31 if strcmpi(options_structure.smooth,'yes') & length(data)==md.numberofelements
    32         data=elementstogrids(md,data);
    33 end
     23%process data and model
     24[x_m y_m z_m elements_m is2d]=processmesh(md,options_structure);
     25[data isongrid]=processdata(md,data,options_structure);
    3426
    35 %layer projection?
    36 if ~isnan(options_structure.layer) & options_structure.layer>=1,
    37         data=project2d(md,data,options_structure.layer); %project onto 2d mesh
    38         %we modify the mesh temporarily to a 2d mesh from which the 3d mesh was extruded.
    39         md.x=md.x2d;
    40         md.y=md.y2d;
    41         md.z=md.z2d;
    42         md.elements=md.elements2d;
    43         md.elements_type=md.elements_type2d;
    44         md.type='2d';
    45 end
    46 
    47 %units
    48 if ~isnan(options_structure.unitmultiplier),
    49         md.x=md.x*options_structure.unitmultiplier;
    50         md.y=md.y*options_structure.unitmultiplier;
    51         md.z=md.z*options_structure.unitmultiplier;
     27%replug x and y onto model so that SectionValue treats the problem correctly
     28if ~isnan(options_structure.layer)
     29        md.x=md.x2d; md.y=md.y2d; md.elements=md.elements2d; md.type='2d';
    5230end
    5331
    … …  
    6139
    6240%Compute section value
    63 [elements,x,y,z,s,data]=SectionValues(md,data,options_structure.sectionvalue,resolution);
     41[elements,x,y,z,s,data_s]=SectionValues(md,data,options_structure.sectionvalue,resolution);
    6442
    65 if strcmpi(md.type,'2d')
     43if is2d
    6644        %plot section value
    6745        subplot(width,width,index1)
    68         plot(s,data)
     46        plot(s,data_s)
     47
     48        %Show Section if requested by user
     49        if ~isnan(options_structure.showsection)
     50
     51                %compute number of labels
     52                numlabels=min(options_structure.showsection,length(s));
     53                shift=fix(length(s)/numlabels);
     54
     55                %plot labels on current graph
     56                hold on
     57                text(s(1),data_s(1),'1','backgroundcolor',[0.8 0.9 0.8])
     58                for i=2:numlabels-1
     59                        text(s(1+(i-1)*shift),data_s(1+(i-1)*shift),num2str(i),'backgroundcolor',[0.8 0.9 0.8])
     60                end
     61                text(s(end),data_s(end),'end','backgroundcolor',[0.8 0.9 0.8])
     62
     63                %plot section only with labels
     64                subplot(width,width,index2)
     65                plot_unit(x_m,y_m,z_m,elements_m,data,isongrid,is2d,options_structure)
     66                hold on
     67                text(x(1),y(1),'1','backgroundcolor',[0.8 0.9 0.8])
     68                for i=2:numlabels-1
     69                        text(x(1+(i-1)*shift),y(1+(i-1)*shift),num2str(i),'backgroundcolor',[0.8 0.9 0.8])
     70                end
     71                text(x(end),y(end),'end','backgroundcolor',[0.8 0.9 0.8])
     72                plot(x,y,'-r')
     73                axis([min(md.x)-1 max(md.x)+1 min(md.y)-1 max(md.y)+1])
     74                view(2)
     75        end
    6976else
    7077        %plot section value
    … …  
    7380                subplot(width,width,index1)
    7481                A=elements(:,1); B=elements(:,2); C=elements(:,3);  D=elements(:,4);
    75                 patch( 'Faces', [A B C D], 'Vertices', [s z zeros(length(s),1)],'FaceVertexCData',data,'FaceColor','interp','EdgeColor','none');
     82                patch( 'Faces', [A B C D], 'Vertices', [s z zeros(length(s),1)],'FaceVertexCData',data_s,'FaceColor','interp','EdgeColor','none');
     83
     84                %Show Section if requested by user
     85                if ~isnan(options_structure.showsection)
     86
     87                        %compute number of labels
     88                        numlabels=min(options_structure.showsection,length(s));
     89                        shift=fix(length(s)/numlabels);
     90
     91                        %plot labels on current graph
     92                        hold on
     93                        text(s(1),z(1),'1','backgroundcolor',[0.8 0.9 0.8])
     94                        for i=2:numlabels-1
     95                                text(s(1+(i-1)*shift),z(1+(i-1)*shift),num2str(i),'backgroundcolor',[0.8 0.9 0.8])
     96                        end
     97                        text(s(end),z(end),'end','backgroundcolor',[0.8 0.9 0.8])
     98
     99                        %plot section only with labels
     100                        subplot(width,width,index2)
     101                        plot_unit(x_m,y_m,z_m,elements_m,data,isongrid,is2d,options_structure)
     102                        hold on
     103                        text(x(1),y(1),'1','backgroundcolor',[0.8 0.9 0.8])
     104                        for i=2:numlabels-1
     105                                text(x(1+(i-1)*shift),y(1+(i-1)*shift),num2str(i),'backgroundcolor',[0.8 0.9 0.8])
     106                        end
     107                        text(x(end),y(end),'end','backgroundcolor',[0.8 0.9 0.8])
     108                        plot(x,y,'-r')
     109                        axis([min(md.x)-1 max(md.x)+1 min(md.y)-1 max(md.y)+1])
     110                        view(2)
     111                end
    76112        else
    77113                subplot(width,width,index1)
    78114                A=elements(:,1); B=elements(:,2); C=elements(:,3);  D=elements(:,4);
    79                 patch( 'Faces', [A B C D], 'Vertices', [x y z],'FaceVertexCData',data,'FaceColor','interp','EdgeColor','none');
     115                patch( 'Faces', [A B C D], 'Vertices', [x y z],'FaceVertexCData',data_s,'FaceColor','interp','EdgeColor','none');
    80116                view(3)
     117
     118                %Show Section if requested by user
     119                if ~isnan(options_structure.showsection)
     120
     121                        %compute number of labels
     122                        numlabels=min(options_structure.showsection,length(s));
     123                        shift=fix(length(x)/numlabels);
     124
     125                        %plot labels on current graph
     126                        hold on
     127                        text(x(1),y(1),z(1),'1','backgroundcolor',[0.8 0.9 0.8])
     128                        for i=2:numlabels-1
     129                                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])
     130                        end
     131                        text(x(end),y(end),z(end),'end','backgroundcolor',[0.8 0.9 0.8])
     132
     133                        %plot section only with labels
     134                        subplot(width,width,index2)
     135                        plot_unit(x_m,y_m,z_m,elements_m,data,isongrid,is2d,options_structure)
     136                        hold on
     137                        text(x(1),y(1),'1','backgroundcolor',[0.8 0.9 0.8])
     138                        for i=2:numlabels-1
     139                                text(x(1+(i-1)*shift),y(1+(i-1)*shift),num2str(i),'backgroundcolor',[0.8 0.9 0.8])
     140                        end
     141                        text(x(end),y(end),'end','backgroundcolor',[0.8 0.9 0.8])
     142                        plot(x,y,'-r')
     143                        axis([min(md.x)-1 max(md.x)+1 min(md.y)-1 max(md.y)+1])
     144                        view(2)
     145                end
    81146        end
    82147end
    … …  
    96161end
    97162applyoptions(md,[],options_structure);
    98 
    99 %plot section if requested by user
    100 if strcmpi(options_structure.showsection,'yes')
    101         subplot(width,width,index2)
    102         hold on
    103         text(x(1),y(1),'1','backgroundcolor',[0.8 0.9 0.8])
    104         text(x(end),y(end),'end','backgroundcolor',[0.8 0.9 0.8])
    105         plot(x,y)
    106         axis([min(md.x) max(md.x) min(md.y) max(md.y)])
    107         view(2)
    108 end

  • issm/trunk/src/m/classes/public/plot/plot_tensor_components.m

    r1 r27  
    55%      plot_tensor_components(md,options_structure,width,i,tensor,type,plot_option);
    66%
    7 %   See also: PLOTMODEL
     7%   See also: PLOTMODEL, PLOT_UNIT, PLOT_DEALER
    88
    9         %Compute the indexes of the components plots
    10         upperplots=fix((i-1)/width);
    11         if upperplots==0, leftplots=i-1; else leftplots=i-width*upperplots-1; end
    12         if strcmpi(md.type,'2d')%3 components -> 3 indexes
    13                 index1=4*width*upperplots+2*leftplots+1;
    14                 index2=index1+1;
    15                 index3=index1+width*2;
    16         elseif strcmpi(md.type,'3d')%6 components -> 6 indexes
    17                 index1=3*3*width*upperplots+3*leftplots+1;
    18                 index2=index1+1;
    19                 index3=index1+2;
    20                 index4=index1+width*3;
    21                 index5=index4+1;
    22                 index6=index4+2;
    23         end
    24 
    25         %smoothing?
    26         if strcmpi(options_structure.smooth,'yes') & length(tensor.xx)==md.numberofelements
    27                 tensor.xx=elementstogrids(md,tensor.xx);
    28                 tensor.yy=elementstogrids(md,tensor.yy);
    29                 tensor.xy=elementstogrids(md,tensor.xy);
    30                 if  strcmpi(md.type,'3d')
    31                 tensor.zz=elementstogrids(md,tensor.zz);
    32                 tensor.yz=elementstogrids(md,tensor.yz);
    33                 tensor.xz=elementstogrids(md,tensor.xz);
    34                 end
    35         end
    36 
    37                        
    38         %layer projection?
    39         if ~isnan(options_structure.layer) & options_structure.layer>=1,
    40                 tensor.xx=project2d(md,tensor.xx,options_structure.layer); %project onto 2d mesh
    41                 tensor.yy=project2d(md,tensor.yy,options_structure.layer); %project onto 2d mesh
    42                 tensor.zz=project2d(md,tensor.zz,options_structure.layer); %project onto 2d mesh
    43                 tensor.xy=project2d(md,tensor.xy,options_structure.layer); %project onto 2d mesh
    44                 tensor.xz=project2d(md,tensor.xz,options_structure.layer); %project onto 2d mesh
    45                 tensor.yz=project2d(md,tensor.yz,options_structure.layer); %project onto 2d mesh
    46                 %we modify the mesh temporarily to a 2d mesh from which the 3d mesh was extruded.
    47                 md.x=md.x2d;
    48                 md.y=md.y2d;
    49                 md.z=md.z2d;
    50                 md.elements=md.elements2d;
    51                 md.elements_type=md.elements_type2d;
    52         end
    53 
    54         %units
    55         if ~isnan(options_structure.unitmultiplier),
    56                 md.x=md.x*options_structure.unitmultiplier;
    57                 md.y=md.y*options_structure.unitmultiplier;
    58                 md.z=md.z*options_structure.unitmultiplier;
    59         end
    60 
    61         if length(tensor.xx)==length(md.elements),
    62                
    63                 if ~isnan(options_structure.noiceshelf) & options_structure.noiceshelf,
    64                         tensor.xx(find(md.elementoniceshelf))=NaN;
    65                         tensor.yy(find(md.elementoniceshelf))=NaN;
    66                         tensor.xy(find(md.elementoniceshelf))=NaN;
    67                         if  strcmpi(md.type,'3d')
    68                                 tensor.zz(find(md.elementoniceshelf))=NaN;
    69                                 tensor.yz(find(md.elementoniceshelf))=NaN;
    70                                 tensor.xz(find(md.elementoniceshelf))=NaN;
    71                         end
    72                 end
    73                 if ~isnan(options_structure.noicesheet) & options_structure.noicesheet,
    74                         tensor.xx(find(~md.elementoniceshelf))=NaN;
    75                         tensor.yy(find(~md.elementoniceshelf))=NaN;
    76                         tensor.xy(find(~md.elementoniceshelf))=NaN;
    77                         if  strcmpi(md.type,'3d')
    78                                 tensor.zz(find(~md.elementoniceshelf))=NaN;
    79                                 tensor.yz(find(~md.elementoniceshelf))=NaN;
    80                                 tensor.xz(find(~md.elementoniceshelf))=NaN;
    81                         end
    82                 end
    83 
    84                 if (strcmpi(md.type,'2d')),
    85                         A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3);
    86 
    87                         subplot(2*width,2*width,index1)
    88                         patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.xx,'FaceColor','flat','EdgeColor','black');
    89                         Apply_options_tensor(options_structure,type,'xx')
    90 
    91                         subplot(2*width,2*width,index2)
    92                         patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.yy,'FaceColor','flat','EdgeColor','black');
    93                         Apply_options_tensor(options_structure,type,'yy')
    94 
    95                         subplot(2*width,2*width,index3)
    96                         patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.xy,'FaceColor','flat','EdgeColor','black');
    97                         Apply_options_tensor(options_structure,type,'xy')
    98                 else
    99                         if options_structure.layer>=1,
    100                                 A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3);
    101 
    102                                 subplot(3*width,3*width,index1)
    103                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData',tensor.xx ,'FaceColor','flat','EdgeColor','black');
    104                                 Apply_options_tensor(options_structure,type,'xx')
    105 
    106                                 subplot(3*width,3*width,index2)
    107                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData',tensor.yy ,'FaceColor','flat','EdgeColor','black');
    108                                 Apply_options_tensor(options_structure,type,'yy')
    109 
    110                                 subplot(3*width,3*width,index3)
    111                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData',tensor.zz ,'FaceColor','flat','EdgeColor','black');
    112                                 Apply_options_tensor(options_structure,type,'zz')
    113 
    114                                 subplot(3*width,3*width,index4)
    115                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData',tensor.xy ,'FaceColor','flat','EdgeColor','black');
    116                                 Apply_options_tensor(options_structure,type,'xy')
    117 
    118                                 subplot(3*width,3*width,index5)
    119                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData',tensor.xz ,'FaceColor','flat','EdgeColor','black');
    120                                 Apply_options_tensor(options_structure,type,'xz')
    121 
    122                                 subplot(3*width,3*width,index6)
    123                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData',tensor.yz ,'FaceColor','flat','EdgeColor','black');
    124                                 Apply_options_tensor(options_structure,type,'yz')
    125 
    126                         else
    127                                 A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3); D=md.elements(:,4); E=md.elements(:,5); F=md.elements(:,6);
    128 
    129                                 subplot(3*width,3*width,index1)
    130                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.xx,'FaceColor','flat','EdgeColor','black');
    131                                 patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'CData', tensor.xx,'FaceColor','flat','EdgeColor','black');
    132                                 patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'CData', tensor.xx,'FaceColor','flat','EdgeColor','black');
    133                                 patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'CData', tensor.xx,'FaceColor','flat','EdgeColor','black');
    134                                 patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'CData', tensor.xx,'FaceColor','flat','EdgeColor','black');
    135                                 Apply_options_tensor(options_structure,type,'xx')
    136 
    137                                 subplot(3*width,3*width,index2)
    138                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.yy,'FaceColor','flat','EdgeColor','black');
    139                                 patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'CData', tensor.yy,'FaceColor','flat','EdgeColor','black');
    140                                 patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'CData', tensor.yy,'FaceColor','flat','EdgeColor','black');
    141                                 patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'CData', tensor.yy,'FaceColor','flat','EdgeColor','black');
    142                                 patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'CData', tensor.yy,'FaceColor','flat','EdgeColor','black');
    143                                 Apply_options_tensor(options_structure,type,'yy')
    144 
    145                                 subplot(3*width,3*width,index3)
    146                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.zz,'FaceColor','flat','EdgeColor','black');
    147                                 patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'CData', tensor.zz,'FaceColor','flat','EdgeColor','black');
    148                                 patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'CData', tensor.zz,'FaceColor','flat','EdgeColor','black');
    149                                 patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'CData', tensor.zz,'FaceColor','flat','EdgeColor','black');
    150                                 patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'CData', tensor.zz,'FaceColor','flat','EdgeColor','black');
    151                                 Apply_options_tensor(options_structure,type,'zz')
    152 
    153                                 subplot(3*width,3*width,index4)
    154                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.xy,'FaceColor','flat','EdgeColor','black');
    155                                 patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'CData', tensor.xy,'FaceColor','flat','EdgeColor','black');
    156                                 patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'CData', tensor.xy,'FaceColor','flat','EdgeColor','black');
    157                                 patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'CData', tensor.xy,'FaceColor','flat','EdgeColor','black');
    158                                 patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'CData', tensor.xy,'FaceColor','flat','EdgeColor','black');
    159                                 Apply_options_tensor(options_structure,type,'xy')
    160 
    161                                 subplot(3*width,3*width,index5)
    162                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.xz,'FaceColor','flat','EdgeColor','black');
    163                                 patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'CData', tensor.xz,'FaceColor','flat','EdgeColor','black');
    164                                 patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'CData', tensor.xz,'FaceColor','flat','EdgeColor','black');
    165                                 patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'CData', tensor.xz,'FaceColor','flat','EdgeColor','black');
    166                                 patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'CData', tensor.xz,'FaceColor','flat','EdgeColor','black');
    167                                 Apply_options_tensor(options_structure,type,'xz')
    168 
    169                                 subplot(3*width,3*width,index6)
    170                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.yz,'FaceColor','flat','EdgeColor','black');
    171                                 patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'CData', tensor.yz,'FaceColor','flat','EdgeColor','black');
    172                                 patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'CData', tensor.yz,'FaceColor','flat','EdgeColor','black');
    173                                 patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'CData', tensor.yz,'FaceColor','flat','EdgeColor','black');
    174                                 patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'CData', tensor.yz,'FaceColor','flat','EdgeColor','black');
    175                                 Apply_options_tensor(options_structure,type,'yz')
    176                         end
    177                 end
    178         elseif length(tensor.xx)==md.numberofgrids |  length(tensor.xx)==md.numberofgrids2d,
    179                 if ~isnan(options_structure.noiceshelf) & options_structure.noiceshelf,
    180                         pos=find(md.gridoniceshelf);
    181                         tensor.xx(pos)=NaN;
    182                         tensor.yy(pos)=NaN;
    183                         tensor.xy(pos)=NaN;
    184                         if strcmpi(md.type,'3d')
    185                                 tensor.zz(pos)=NaN;
    186                                 tensor.xz(pos)=NaN;
    187                                 tensor.yz(pos)=NaN;
    188                         end
    189                 end
    190                 if ~isnan(options_structure.noiceshelf) & options_structure.noicesheet,
    191                         pos=find(md.gridonicesheet);
    192                         tensor.xx(pos)=NaN;
    193                         tensor.yy(pos)=NaN;
    194                         tensor.xy(pos)=NaN;
    195                         if strcmpi(md.type,'3d')
    196                                 tensor.zz(pos)=NaN;
    197                                 tensor.xz(pos)=NaN;
    198                                 tensor.yz(pos)=NaN;
    199                         end
    200                 end
    201 
    202                 if strcmpi(md.type,'2d'),
    203                         A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3);
    204 
    205                         subplot(2*width,2*width,index1)
    206                         patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.xx ,'FaceColor','flat','EdgeColor','none');
    207                         Apply_options_tensor(options_structure,type,'xx')
    208 
    209                         subplot(2*width,2*width,index2)
    210                         patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.yy ,'FaceColor','flat','EdgeColor','none');
    211                         Apply_options_tensor(options_structure,type,'yy')
    212 
    213                         subplot(2*width,2*width,index3)
    214                         patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.xy ,'FaceColor','flat','EdgeColor','none');
    215                         Apply_options_tensor(options_structure,type,'xy')
    216 
    217                 else
    218                         if options_structure.layer>=1,
    219                                 A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3);
    220                                 subplot(3*width,3*width,index1)
    221                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.xx ,'FaceColor','flat','EdgeColor','none');
    222                                 Apply_options_tensor(options_structure,type,'xx')
    223 
    224                                 subplot(3*width,3*width,index2)
    225                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.yy ,'FaceColor','flat','EdgeColor','none');
    226                                 Apply_options_tensor(options_structure,type,'yy')
    227 
    228                                 subplot(3*width,3*width,index3)
    229                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.zz ,'FaceColor','flat','EdgeColor','none');
    230                                 Apply_options_tensor(options_structure,type,'zz')
    231 
    232                                 subplot(3*width,3*width,index4)
    233                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.xy ,'FaceColor','flat','EdgeColor','none');
    234                                 Apply_options_tensor(options_structure,type,'xy')
    235 
    236                                 subplot(3*width,3*width,index5)
    237                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.xz ,'FaceColor','flat','EdgeColor','none');
    238                                 Apply_options_tensor(options_structure,type,'xz')
    239 
    240                                 subplot(3*width,3*width,index6)
    241                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.yz ,'FaceColor','flat','EdgeColor','none');
    242                                 Apply_options_tensor(options_structure,type,'yz')
    243 
    244                         else
    245                                 A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3); D=md.elements(:,4); E=md.elements(:,5); F=md.elements(:,6);
    246 
    247                                 subplot(3*width,3*width,index1)
    248                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xx,'FaceColor','flat','EdgeColor','none');
    249                                 patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xx,'FaceColor','flat','EdgeColor','none');
    250                                 patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xx,'FaceColor','flat','EdgeColor','none');
    251                                 patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xx,'FaceColor','flat','EdgeColor','none');
    252                                 patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xx,'FaceColor','flat','EdgeColor','none');
    253                                 Apply_options_tensor(options_structure,type,'xx')
    254 
    255                                 subplot(3*width,3*width,index2)
    256                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.yy,'FaceColor','flat','EdgeColor','none');
    257                                 patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.yy,'FaceColor','flat','EdgeColor','none');
    258                                 patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.yy,'FaceColor','flat','EdgeColor','none');
    259                                 patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.yy,'FaceColor','flat','EdgeColor','none');
    260                                 patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.yy,'FaceColor','flat','EdgeColor','none');
    261                                 Apply_options_tensor(options_structure,type,'yy')
    262 
    263                                 subplot(3*width,3*width,index3)
    264                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.zz,'FaceColor','flat','EdgeColor','none');
    265                                 patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.zz,'FaceColor','flat','EdgeColor','none');
    266                                 patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.zz,'FaceColor','flat','EdgeColor','none');
    267                                 patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.zz,'FaceColor','flat','EdgeColor','none');
    268                                 patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.zz,'FaceColor','flat','EdgeColor','none');
    269                                 Apply_options_tensor(options_structure,type,'zz')
    270 
    271                                 subplot(3*width,3*width,index4)
    272                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xy,'FaceColor','flat','EdgeColor','none');
    273                                 patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xy,'FaceColor','flat','EdgeColor','none');
    274                                 patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xy,'FaceColor','flat','EdgeColor','none');
    275                                 patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xy,'FaceColor','flat','EdgeColor','none');
    276                                 patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xy,'FaceColor','flat','EdgeColor','none');
    277                                 Apply_options_tensor(options_structure,type,'xy')
    278 
    279                                 subplot(3*width,3*width,index5)
    280                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xz,'FaceColor','flat','EdgeColor','none');
    281                                 patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xz,'FaceColor','flat','EdgeColor','none');
    282                                 patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xz,'FaceColor','flat','EdgeColor','none');
    283                                 patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xz,'FaceColor','flat','EdgeColor','none');
    284                                 patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.xz,'FaceColor','flat','EdgeColor','none');
    285                                 Apply_options_tensor(options_structure,type,'xz')
    286 
    287                                 subplot(3*width,3*width,index6)
    288                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.yz,'FaceColor','flat','EdgeColor','none');
    289                                 patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.yz,'FaceColor','flat','EdgeColor','none');
    290                                 patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.yz,'FaceColor','flat','EdgeColor','none');
    291                                 patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.yz,'FaceColor','flat','EdgeColor','none');
    292                                 patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.yz,'FaceColor','flat','EdgeColor','none');
    293                                 Apply_options_tensor(options_structure,type,'yz')
    294                         end
    295                 end
    296         end
     9%Compute the indexes of the components plots
     10upperplots=fix((i-1)/width);
     11if upperplots==0, leftplots=i-1; else leftplots=i-width*upperplots-1; end
     12if strcmpi(md.type,'2d')%3 components -> 3 indexes
     13        index1=4*width*upperplots+2*leftplots+1;
     14        index2=index1+1;
     15        index3=index1+width*2;
     16elseif strcmpi(md.type,'3d')%6 components -> 6 indexes
     17        index1=3*3*width*upperplots+3*leftplots+1;
     18        index2=index1+1;
     19        index3=index1+2;
     20        index4=index1+width*3;
     21        index5=index4+1;
     22        index6=index4+2;
    29723end
    29824
    299 function Apply_options_tensor(options_structure,type,component)
    300 %apply options
     25%process data and model
     26[x y z elements is2d]=processmesh(md,options_structure);
     27[tensor.xx isongrid]=processdata(md,tensor.xx,options_structure);
     28[tensor.yy isongrid]=processdata(md,tensor.yy,options_structure);
     29[tensor.xy isongrid]=processdata(md,tensor.xy,options_structure);
     30if  strcmpi(md.type,'3d')
     31        [tensor.xz isongrid]=processdata(md,tensor.xz,options_structure);
     32        [tensor.yz isongrid]=processdata(md,tensor.yz,options_structure);
     33        [tensor.zz isongrid]=processdata(md,tensor.zz,options_structure);
     34end
     35
     36if (strcmpi(md.type,'2d')),
     37        subplot(2*width,2*width,index1),
     38        plot_unit(x,y,z,elements,tensor.xx,isongrid,is2d,options_structure)
     39        Apply_options_tensor(md,options_structure,type,'xx')
     40        subplot(2*width,2*width,index2),
     41        plot_unit(x,y,z,elements,tensor.yy,isongrid,is2d,options_structure)
     42        Apply_options_tensor(md,options_structure,type,'yy')
     43        subplot(2*width,2*width,index3),
     44        plot_unit(x,y,z,elements,tensor.xy,isongrid,is2d,options_structure)
     45        Apply_options_tensor(md,options_structure,type,'xy')
     46else
     47        subplot(3*width,3*width,index1),
     48        plot_unit(x,y,z,elements,tensor.xx,isongrid,is2d,options_structure)
     49        Apply_options_tensor(md,options_structure,type,'xx')
     50        subplot(3*width,3*width,index2),
     51        plot_unit(x,y,z,elements,tensor.yy,isongrid,is2d,options_structure)
     52        Apply_options_tensor(md,options_structure,type,'yy')
     53        subplot(3*width,3*width,index3),
     54        plot_unit(x,y,z,elements,tensor.zz,isongrid,is2d,options_structure)
     55        Apply_options_tensor(md,options_structure,type,'zz')
     56        subplot(3*width,3*width,index4),
     57        plot_unit(x,y,z,elements,tensor.xy,isongrid,is2d,options_structure)
     58        Apply_options_tensor(md,options_structure,type,'xy')
     59        subplot(3*width,3*width,index5),
     60        plot_unit(x,y,z,elements,tensor.xz,isongrid,is2d,options_structure)
     61        Apply_options_tensor(md,options_structure,type,'xz')
     62        subplot(3*width,3*width,index6),
     63        plot_unit(x,y,z,elements,tensor.yz,isongrid,is2d,options_structure)
     64        Apply_options_tensor(md,options_structure,type,'yz')
     65end
     66end
     67
     68function Apply_options_tensor(md,options_structure,type,component)
     69        %apply options
    30170        if isnan(options_structure.title)
    30271                if ismember('_',type) %user plotet stress_tensor

  • issm/trunk/src/m/classes/public/plot/plot_tensor_principal.m

    r1 r27  
    55%      plot_tensor_principal(md,options_structure,width,i,tensor,type,plot_options);
    66%
    7 %   See also: PLOTMODEL
     7%   See also: PLOTMODEL, PLOT_UNIT, PLOT_DEALER
    88
    9         %Compute the indexes of the components plots
    10         upperplots=fix((i-1)/width);
    11         if upperplots==0, leftplots=i-1; else leftplots=i-width*upperplots-1; end
    12         if strcmpi(md.type,'2d')%3 components -> 3 indexes
    13                 index1=4*width*upperplots+2*leftplots+1;
    14                 index2=index1+1;
    15                 index3=index1+width*2;
    16                 index4=index3+1;
    17                 newwidth=2*width;
    18         elseif strcmpi(md.type,'3d')%6 components -> 6 indexes
    19                 index1=3*3*width*upperplots+3*leftplots+1;
    20                 index2=index1+1;
    21                 index3=index1+2;
    22                 index4=index1+width*3;
    23                 index5=index4+1;
    24                 index6=index4+2;
    25                 newwidth=3*width;
    26         end
    27 
    28         %plot principal axis
    29         type1=[type 'axis1'];
    30         plot_tensor_principalaxis(md,options_structure,newwidth,index1,tensor,type1,plot_options);
    31         type2=[type 'axis2'];
    32         plot_tensor_principalaxis(md,options_structure,newwidth,index2,tensor,type2,plot_options);
    33         if  strcmpi(md.type,'3d')
    34                 type3=[type 'axis3'];
    35                 plot_tensor_principalaxis(md,options_structure,newwidth,index3,tensor,type3,plot_options);
    36         end
    37 
    38         %smoothing?
    39         if strcmpi(options_structure.smooth,'yes') & length(tensor.principalvalue1)==md.numberofelements
    40                 tensor.principalvalue1=elementstogrids(md,tensor.principalvalue1);
    41                 tensor.principalvalue2=elementstogrids(md,tensor.principalvalue2);
    42                 if  strcmpi(md.type,'3d'), tensor.principalvalue3=elementstogrids(md,tensor.principalvalue3); end
    43         end
    44 
    45                        
    46         %layer projection?
    47         if ~isnan(options_structure.layer) & options_structure.layer>=1,
    48                 tensor.principalvalue1=project2d(md,tensor.principalvalue1,options_structure.layer); %project onto 2d mesh
    49                 tensor.principalvalue2=project2d(md,tensor.principalvalue2,options_structure.layer);
    50                 tensor.principalvalue3=project2d(md,tensor.principalvalue3,options_structure.layer);
    51 
    52                 %we modify the mesh temporarily to a 2d mesh from which the 3d mesh was extruded.
    53                 md.x=md.x2d;
    54                 md.y=md.y2d;
    55                 md.z=md.z2d;
    56                 md.elements=md.elements2d;
    57                 md.elements_type=md.elements_type2d;
    58         end
    59 
    60         %units
    61         if ~isnan(options_structure.unitmultiplier),
    62                 md.x=md.x*options_structure.unitmultiplier;
    63                 md.y=md.y*options_structure.unitmultiplier;
    64                 md.z=md.z*options_structure.unitmultiplier;
    65         end
    66 
    67         if length(tensor.principalvalue1)==length(md.elements),
    68                 if ~isnan(options_structure.noiceshelf) & options_structure.noiceshelf,
    69                         tensor.principalvalue1(find(md.elementoniceshelf))=NaN;
    70                         tensor.principalvalue2(find(md.elementoniceshelf))=NaN;
    71                         if  strcmpi(md.type,'3d')
    72                                 tensor.principalvalue3(find(md.elementoniceshelf))=NaN;
    73                         end
    74                 end
    75                 if ~isnan(options_structure.noicesheet) & options_structure.noicesheet,
    76                         tensor.principalvalue1(find(~md.elementoniceshelf))=NaN;
    77                         tensor.principalvalue2(find(~md.elementoniceshelf))=NaN;
    78                         if  strcmpi(md.type,'3d')
    79                                 tensor.principalvalue3(find(~md.elementoniceshelf))=NaN;
    80                         end
    81                 end
    82 
    83                 if (strcmpi(md.type,'2d')),
    84                         A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3);
    85 
    86                         subplot(2*width,2*width,index3)
    87                         patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue1,'FaceColor','flat','EdgeColor','black');
    88                         Apply_options_tensor(options_structure,type,'principal value 1')
    89 
    90                         subplot(2*width,2*width,index4)
    91                         patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue2,'FaceColor','flat','EdgeColor','black');
    92                         Apply_options_tensor(options_structure,type,'principal value 2')
    93                 else
    94                         if options_structure.layer>=1,
    95                                 A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3);
    96 
    97                                 subplot(3*width,3*width,index4)
    98                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue1,'FaceColor','flat','EdgeColor','black');
    99                                 Apply_options_tensor(options_structure,type,'principal value 1')
    100 
    101                                 subplot(3*width,3*width,index5)
    102                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue2,'FaceColor','flat','EdgeColor','black');
    103                                 Apply_options_tensor(options_structure,type,'principal value 2')
    104 
    105                                 subplot(3*width,3*width,index6)
    106                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue3,'FaceColor','flat','EdgeColor','black');
    107                                 Apply_options_tensor(options_structure,type,'principal value 3')
    108 
    109                         else
    110        
    111                                 A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3); D=md.elements(:,4); E=md.elements(:,5); F=md.elements(:,6);
    112 
    113                                 subplot(3*width,3*width,index4)
    114                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue1,'FaceColor','flat','EdgeColor','black');
    115                                 patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue1,'FaceColor','flat','EdgeColor','black');
    116                                 patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue1,'FaceColor','flat','EdgeColor','black');
    117                                 patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue1,'FaceColor','flat','EdgeColor','black');
    118                                 patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue1,'FaceColor','flat','EdgeColor','black');
    119                                 Apply_options_tensor(options_structure,type,'principal value 1')
    120 
    121                                 subplot(3*width,3*width,index5)
    122                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue2,'FaceColor','flat','EdgeColor','black');
    123                                 patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue2,'FaceColor','flat','EdgeColor','black');
    124                                 patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue2,'FaceColor','flat','EdgeColor','black');
    125                                 patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue2,'FaceColor','flat','EdgeColor','black');
    126                                 patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue2,'FaceColor','flat','EdgeColor','black');
    127                                 Apply_options_tensor(options_structure,type,'principal value 2')
    128 
    129                                 subplot(3*width,3*width,index6)
    130                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue3,'FaceColor','flat','EdgeColor','black');
    131                                 patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue3,'FaceColor','flat','EdgeColor','black');
    132                                 patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue3,'FaceColor','flat','EdgeColor','black');
    133                                 patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue3,'FaceColor','flat','EdgeColor','black');
    134                                 patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'CData', tensor.principalvalue3,'FaceColor','flat','EdgeColor','black');
    135                                 Apply_options_tensor(options_structure,type,'principal value 3')
    136                         end
    137                 end
    138 
    139         elseif length(tensor.principalvalue1)==md.numberofgrids |  length(tensor.principalvalue1)==md.numberofgrids2d,
    140                 if ~isnan(options_structure.noiceshelf) & options_structure.noiceshelf,
    141                         pos=find(md.gridoniceshelf);
    142                         tensor.principalvalue1(pos)=NaN;
    143                         tensor.principalvalue2(pos)=NaN;
    144                         if strcmpi(md.type,'3d')
    145                                 tensor.principalvalue3(pos)=NaN;
    146                         end
    147                 end
    148                 if ~isnan(options_structure.noiceshelf) & options_structure.noicesheet,
    149                         pos=find(md.gridonicesheet);
    150                         tensor.principalvalue1(pos)=NaN;
    151                         tensor.principalvalue2(pos)=NaN;
    152                         if strcmpi(md.type,'3d')
    153                                 tensor.principalvalue3(pos)=NaN;
    154                         end
    155                 end
    156 
    157                 if strcmpi(md.type,'2d'),
    158                         A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3);
    159 
    160                         subplot(2*width,2*width,index3)
    161                         patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.principalvalue1 ,'FaceColor','flat','EdgeColor','none');
    162                         Apply_options_tensor(options_structure,type,'principal value 1')
    163 
    164                         subplot(2*width,2*width,index4)
    165                         patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.principalvalue2 ,'FaceColor','flat','EdgeColor','none');
    166                         Apply_options_tensor(options_structure,type,'principal value 2')
    167                 else
    168                         if options_structure.layer>=1,
    169                                 A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3);
    170 
    171                                 subplot(3*width,3*width,index4)
    172                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.principalvalue1 ,'FaceColor','flat','EdgeColor','none');
    173                                 Apply_options_tensor(options_structure,type,'principal value 1')
    174 
    175                                 subplot(3*width,3*width,index5)
    176                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.principalvalue2 ,'FaceColor','flat','EdgeColor','none');
    177                                 Apply_options_tensor(options_structure,type,'principal value 2')
    178 
    179                                 subplot(3*width,3*width,index6)
    180                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData',tensor.principalvalue3 ,'FaceColor','flat','EdgeColor','none');
    181                                 Apply_options_tensor(options_structure,type,'principal value 3')
    182                         else
    183                                 A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3); D=md.elements(:,4); E=md.elements(:,5); F=md.elements(:,6);
    184 
    185                                 subplot(3*width,3*width,index4)
    186                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue1,'FaceColor','flat','EdgeColor','none');
    187                                 patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue1,'FaceColor','flat','EdgeColor','none');
    188                                 patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue1,'FaceColor','flat','EdgeColor','none');
    189                                 patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue1,'FaceColor','flat','EdgeColor','none');
    190                                 patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue1,'FaceColor','flat','EdgeColor','none');
    191                                 Apply_options_tensor(options_structure,type,'principal value 1')
    192 
    193                                 subplot(3*width,3*width,index5)
    194                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue2,'FaceColor','flat','EdgeColor','none');
    195                                 patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue2,'FaceColor','flat','EdgeColor','none');
    196                                 patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue2,'FaceColor','flat','EdgeColor','none');
    197                                 patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue2,'FaceColor','flat','EdgeColor','none');
    198                                 patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue2,'FaceColor','flat','EdgeColor','none');
    199                                 Apply_options_tensor(options_structure,type,'principal value 2')
    200 
    201                                 subplot(3*width,3*width,index6)
    202                                 patch( 'Faces', [A B C], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue3,'FaceColor','flat','EdgeColor','none');
    203                                 patch( 'Faces', [D E F], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue3,'FaceColor','flat','EdgeColor','none');
    204                                 patch( 'Faces', [A B E D], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue3,'FaceColor','flat','EdgeColor','none');
    205                                 patch( 'Faces', [B E F C ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue3,'FaceColor','flat','EdgeColor','none');
    206                                 patch( 'Faces', [C A D F ], 'Vertices', [md.x md.y md.z],'FaceVertexCData', tensor.principalvalue3,'FaceColor','flat','EdgeColor','none');
    207                                 Apply_options_tensor(options_structure,type,'principal value 3')
    208                         end
    209                 end
    210         end
     9%Compute the indexes of the components plots
     10upperplots=fix((i-1)/width);
     11if upperplots==0, leftplots=i-1; else leftplots=i-width*upperplots-1; end
     12if strcmpi(md.type,'2d')%3 components -> 3 indexes
     13        index1=4*width*upperplots+2*leftplots+1;
     14        index2=index1+1;
     15        index3=index1+width*2;
     16        index4=index3+1;
     17        newwidth=2*width;
     18elseif strcmpi(md.type,'3d')%6 components -> 6 indexes
     19        index1=3*3*width*upperplots+3*leftplots+1;
     20        index2=index1+1;
     21        index3=index1+2;
     22        index4=index1+width*3;
     23        index5=index4+1;
     24        index6=index4+2;
     25        newwidth=3*width;
    21126end
    21227
    213 function Apply_options_tensor(options_structure,type,component)
     28%plot principal axis
     29type1=[type 'axis1'];
     30plot_tensor_principalaxis(md,options_structure,newwidth,index1,tensor,type1,plot_options);
     31type2=[type 'axis2'];
     32plot_tensor_principalaxis(md,options_structure,newwidth,index2,tensor,type2,plot_options);
     33if  strcmpi(md.type,'3d')
     34        type3=[type 'axis3'];
     35        plot_tensor_principalaxis(md,options_structure,newwidth,index3,tensor,type3,plot_options);
     36end
     37
     38%now plot principal values
     39[x y z elements is2d]=processmesh(md,options_structure);
     40[tensor.principalvalue1 isongrid]=processdata(md,tensor.principalvalue1,options_structure);
     41[tensor.principalvalue2 isongrid]=processdata(md,tensor.principalvalue2,options_structure);
     42if  strcmpi(md.type,'3d')
     43        [tensor.principalvalue3 isongrid]=processdata(md,tensor.principalvalue3,options_structure);
     44end
     45
     46if (strcmpi(md.type,'2d')),
     47        subplot(2*width,2*width,index3)
     48        plot_unit(x,y,z,elements,tensor.principalvalue1,isongrid,is2d,options_structure)
     49        Apply_options_tensor(md,options_structure,type,'principal value 1')
     50        subplot(2*width,2*width,index4)
     51        plot_unit(x,y,z,elements,tensor.principalvalue2,isongrid,is2d,options_structure)
     52        Apply_options_tensor(md,options_structure,type,'principal value 2')
     53else
     54        subplot(3*width,3*width,index4)
     55        plot_unit(x,y,z,elements,tensor.principalvalue1,isongrid,is2d,options_structure)
     56        Apply_options_tensor(md,options_structure,type,'principal value 1')
     57        subplot(3*width,3*width,index5)
     58        plot_unit(x,y,z,elements,tensor.principalvalue2,isongrid,is2d,options_structure)
     59        Apply_options_tensor(md,options_structure,type,'principal value 2')
     60        subplot(3*width,3*width,index6)
     61        plot_unit(x,y,z,elements,tensor.principalvalue3,isongrid,is2d,options_structure)
     62        Apply_options_tensor(md,options_structure,type,'principal value 3')
     63end
     64end
     65
     66function Apply_options_tensor(md,options_structure,type,component)
    21467%apply options
    21568        if isnan(options_structure.title)

  • issm/trunk/src/m/classes/public/plot/plot_tensor_principalaxis.m

    r1 r27  
    77%   See also: PLOTMODEL
    88
    9 %plot mesh boundaries
     9%prepare subplot
    1010subplot(width,width,i);
    1111
     12%process data and model
     13[x y z elements is2d]=processmesh(md,options_structure);
    1214
    1315if (strcmpi(md.type,'2d')),
    1416        eval(['Vx=tensor.principalaxis' type(end) '(:,1); Vy=tensor.principalaxis' type(end) '(:,2);'])
     17        eval(['value=tensor.principalvalue' type(end) ';']);
     18        [Vx isongrid]=processdata(md,Vx,options_structure);
     19        [Vy isongrid]=processdata(md,Vy,options_structure);
     20        [value isongrid]=processdata(md,value,options_structure);
    1521else
    1622        eval(['Vx=tensor.principalaxis' type(end) '(:,1); Vy=tensor.principalaxis' type(end) '(:,2); Vz=tensor.principalaxis' type(end) '(:,3);'])
     23        [Vx isongrid]=processdata(md,Vx,options_structure);
     24        [Vy isongrid]=processdata(md,Vy,options_structure);
     25        [Vz isongrid]=processdata(md,Vz,options_structure);
     26        [value isongrid]=processdata(md,value,options_structure);
    1727end
    1828
    19 %smoothing?
    20 if strcmpi(options_structure.smooth,'yes') & length(Vx)==md.numberofelements
    21         Vx=elementstogrids(md,Vx);
    22         Vy=elementstogrids(md,Vy);
    23         if (strcmpi(md.type,'3d')),
    24                 Vz=elementstogrids(md,Vz);
    25         end
    26 end
    27                
    28 %layer projection?
    29 if ~isnan(options_structure.layer) & options_structure.layer>=1,
    30         Vx=project2d(md,Vx,options_structure.layer); %project onto 2d mesh
    31         Vy=project2d(md,Vy,options_structure.layer); %project onto 2d mesh
    32         if (strcmpi(md.type,'3d')),
    33                 Vz=project2d(md,Vz,options_structure.layer); %project onto 2d mesh
    34         end
    35         %we modify the mesh temporarily to a 2d mesh from which the 3d mesh was extruded.
    36         md.x=md.x2d;
    37         md.y=md.y2d;
    38         md.z=md.z2d;
    39         md.elements=md.elements2d;
    40         md.elements_type=md.elements_type2d;
    41         md.type='2d';
     29%take the center of each element if ~isongrid
     30if ~isongrid
     31        x=mean(md.x(md.elements'))'; y=mean(md.y(md.elements'))'; z=mean(md.z(md.elements'))';
    4232end
    4333
    44 %units
    45 if ~isnan(options_structure.unitmultiplier),
    46         md.x=md.x*options_structure.unitmultiplier;
    47         md.y=md.y*options_structure.unitmultiplier;
    48         md.z=md.z*options_structure.unitmultiplier;
    49 end
    50 if length(Vx)==length(md.elements),
    51        
    52         if ~isnan(options_structure.noiceshelf) & options_structure.noiceshelf,
    53                 Vx(find(md.elementoniceshelf))=NaN;
    54                 Vy(find(md.elementoniceshelf))=NaN;
    55                 Vz(find(md.elementoniceshelf))=NaN;
    56         end
    57         if ~isnan(options_structure.noicesheet) & options_structure.noicesheet,
    58                 Vx(find(~md.elementoniceshelf))=NaN;
    59                 Vy(find(~md.elementoniceshelf))=NaN;
    60                 Vz(find(~md.elementoniceshelf))=NaN;
     34%plot quivers
     35if strcmpi(md.type,'2d'),
     36
     37        %density
     38        if ~isnan(options_structure.density)
     39                x=x(1:options_structure.density:end);
     40                y=y(1:options_structure.density:end);
     41                Vx=Vx(1:options_structure.density:end);
     42                Vy=Vy(1:options_structure.density:end);
     43                value=value(1:options_structure.density:end);
    6144        end
    6245
    63         if (strcmpi(md.type,'2d')),
    64                 x=mean(md.x(md.elements'))'; y=mean(md.y(md.elements'))';
    65                 A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3);
    66                 quiver(x,y,Vx,Vy)
     46        %scaling:
     47        delta=((min(x)-max(x))^2+(min(y)-max(y))^2)/numel(x);
     48        scale=0.5/max(sqrt((Vx.^2+Vy.^2)/delta));
     49        Vx=scale*Vx; Vy=scale*Vy;
    6750
    68         else
    69                 x=mean(md.x(md.elements'))'; y=mean(md.y(md.elements'))'; z=mean(md.z(md.elements'))';
    70                 A=md.elements(:,1); B=md.elements(:,2); C=md.elements(:,3); D=md.elements(:,4); E=md.elements(:,5); F=md.elements(:,6);
    71                 quiver3(x,y,z,Vx,Vy,Vz)
    72         end
    73 elseif length(Vx)==md.numberofgrids |  length(Vx)==md.numberofgrids2d,
    74         if ~isnan(options_structure.noiceshelf) & options_structure.noiceshelf,
    75                 pos=find(md.gridoniceshelf);
    76                 Vx(pos)=NaN;
    77                 Vy(pos)=NaN;
    78                 Vz(pos)=NaN;
    79         end
    80         if ~isnan(options_structure.noiceshelf) & options_structure.noicesheet,
    81                 pos=find(md.gridonicesheet);
    82                 Vx(pos)=NaN;
    83                 Vy(pos)=NaN;
    84                 Vz(pos)=NaN;
     51        pos=find(value>=0);
     52        q1=quiver(x(pos),y(pos),Vx(pos),Vy(pos),'Color','r','ShowArrowHead','off','AutoScale','off');
     53        hold on
     54        pos=find(value<0);
     55        q2=quiver(x(pos),y(pos),Vx(pos),Vy(pos),'Color','b','ShowArrowHead','off','AutoScale','off');
     56
     57else
     58        %density
     59        if ~isnan(options_structure.density)
     60                x=x(1:options_structure.density:end);
     61                y=y(1:options_structure.density:end);
     62                z=z(1:options_structure.density:end);
     63                Vx=Vx(1:options_structure.density:end);
     64                Vy=Vy(1:options_structure.density:end);
     65                Vz=Vz(1:options_structure.density:end);
     66                value=value(1:options_structure.density:end);
    8567        end
    8668
    87         if strcmpi(md.type,'2d'),
    88                 quiver(md.x,md.y,Vx,Vy)
    89         else
    90                 quiver3(md.x,md.y,md.z,Vx,Vy,Vz)
    91         end
     69        %scaling:
     70        delta=((min(x)-max(x))^2+(min(y)-max(y))^2)/numel(x);
     71        scale=0.5/max(sqrt((Vx.^2+Vy.^2)/delta));
     72        Vx=scale*Vx; Vy=scale*Vy; Vz=scale*Vz;
     73
     74        pos=find(value>=0);
     75        q1=quiver3(x(pos),y(pos),z(pos),Vx(pos),Vy(pos),Vz(pos),'Color','r','ShowArrowHead','off','AutoScale','off');
     76        hold on
     77        pos=find(value<0);
     78        q2=quiver3(x(pos),y(pos),z(pos),Vx(pos),Vy(pos),Vz(pos),'Color','b','ShowArrowHead','off','AutoScale','off');
     79end
     80
     81%legend
     82if strcmpi(type(1:6),'strain')
     83        legend([q1 q2],'extension','compression')
     84elseif strcmpi(type(1:6),'stress')
     85        legend([q1 q2],'compression','traction')
    9286end
    9387

  • issm/trunk/src/m/classes/public/plot/plot_transient_movie.m

    r1 r27  
    44%      plot_transient_movie(md,options_structure,width,i);
    55%
    6 %   See also: PLOTMODEL
     6%   See also: PLOTMODEL, PLOT_UNIT, PLOT_DEALER
    77
    8         %plot mesh boundaries
     8        %prepare subplot
    99        subplot(width,width,i);
    10 
    11         %first load x,y, etc ... to speed up plot
    12         x=md.x;
    13         x2d=md.x2d;
    14         y=md.y;
    15         y2d=md.y2d;
    16         z=md.z;
    17         z2d=md.z2d;
    18         elements2d=md.elements2d;
    19         elements=md.elements;
    20         elements_type2d=md.elements_type2d;
    21 
    22         %units
    23         if ~isnan(options_structure.unitmultiplier),
    24                 md.x=md.x*options_structure.unitmultiplier;
    25                 md.y=md.y*options_structure.unitmultiplier;
    26                 md.z=md.z*options_structure.unitmultiplier;
    27         end
    28 
    29         %edgecolor?
    30         if ~isnan(options_structure.edgecolor),
    31                 edgecolor=options_structure.edgecolor;
    32         else
    33                 edgecolor='none';
    34         end
    3510
    3611        if strcmpi(md.type,'2d')
    3712                choice=input('Which field do you want to plot? (vel/vx/vy/thickness/surface/bed)','s');
    38 
    3913                if ~strcmp(choice,'vel') & ~strcmp(choice,'vx') & ~strcmp(choice,'vy') & ~strcmp(choice,'thickness') & ~strcmp(choice,'bed') & ~strcmp(choice,'surface')
    4014                        disp('plot_transient_movie error message: input not supported yet, exiting...')
    4115                        return
    4216                end
    43                 if ~isnan(options_structure.noiceshelf) & options_structure.noiceshelf,
    44                         pos=find(md.gridoniceshelf);
    45                         data(pos)=NaN;
    46                 end
    47                 if ~isnan(options_structure.noiceshelf) & options_structure.noicesheet,
    48                         pos=find(md.gridonicesheet);
    49                         data(pos)=NaN;
    50                 end
    51                 for i=1:length(md.transient_results)
    52                         eval(['data=md.transient_results(' num2str(i) ').' num2str(choice) ';']);
    53                         titlestring=[choice ' at time ' num2str(md.transient_results(i).time) ' year'];
    54                         A=elements(:,1); B=elements(:,2); C=elements(:,3);
    55                         patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor);
    56                         apply_options_movie(options_structure,titlestring);
    57                         pause(0.5)
    58                 end
    59 
    6017        else
    6118                choice=input('Which field do you want to plot? (vel/vx/vy/vz/thickness/surface/bed/temperature)','s');
    62 
    6319                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')
    6420                        disp('plot_transient_movie error message: input not supported yet, exiting...')
    6521                        return
    6622                end
    67                 if ~isnan(options_structure.noiceshelf) & options_structure.noiceshelf,
    68                         pos=find(md.gridoniceshelf);
    69                         data(pos)=NaN;
    70                 end
    71                 if ~isnan(options_structure.noiceshelf) & options_structure.noicesheet,
    72                         pos=find(md.gridonicesheet);
    73                         data(pos)=NaN;
    74                 end
    75                 for i=1:length(md.transient_results)
    76                         eval(['data=md.transient_results(' num2str(i) ').' num2str(choice) ';']);
    77                         titlestring=[choice 'at time ' num2str(md.transient_results(i).time) ' a'];
    78                         A=elements(:,1); B=elements(:,2); C=elements(:,3); D=elements(:,4); E=elements(:,5); F=elements(:,6);
    79                         patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor);
    80                         patch( 'Faces', [D E F], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor);
    81                         patch( 'Faces', [A B E D], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor);
    82                         patch( 'Faces', [B E F C ], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor);
    83                         patch( 'Faces', [C A D F ], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor);
    84                         apply_options_movie(options_structure,titlestring);
    85                         pause(0.5)
    86                 end
    8723        end
    8824
     25        %process model
     26        [x y z elements is2d]=processmesh(md,options_structure);
     27
     28        %loop over the time steps
     29        for i=1:length(md.transient_results)
     30                eval(['data=md.transient_results(' num2str(i) ').' num2str(choice) ';']);
     31
     32                %process data
     33                [data isongrid]=processdata(md,data,options_structure);
     34                titlestring=[choice ' at time ' num2str(md.transient_results(i).time) ' year'];
     35                plot_unit(x,y,z,elements,data,isongrid,is2d,options_structure)
     36                apply_options_movie(md,options_structure,titlestring);
     37                pause(0.5)
     38        end
    8939end %function
    9040
    91 function apply_options_movie(options_structure,titlestring)
     41function apply_options_movie(md,options_structure,titlestring)
    9242        %apply options
    9343        if isnan(options_structure.title)

  • issm/trunk/src/m/classes/public/plot/plot_unit.m

    r1 r27  
    1 function plot_unit(md,optionstring,width,i);
    2 %PLOT_UNIT - unit plot called by plotmodel
     1function plot_unit(x,y,z,elements,data,isongrid,is2d,options_structure)
     2%PLOT_UNIT - unit plot, display data
    33%
    44%   Usage:
    5 %      plot_unit(md,optionstring,width,i);
     5%      plot_unit(x,y,z,elements,data,isongrid,is2d,options_structure);
    66%
    7 %   See also: PLOTMODEL
    8 
    9 %parse options and get a structure of options.
    10 options_structure=parse_options(md,optionstring);
    11 
    12 %figure out if this is a special plot
    13 data=findarg(optionstring,'data');data=data.value;
    14 
    15 %Figure out if this is a semi-transparent plot.
    16 if ~isnan(options_structure.overlay),
    17         plot_overlay(md,data,options_structure,width,i);
    18         return;
    19 end
    20 
    21 if ischar(data),
    22         switch data,
    23         case 'boundaries',
    24                 plot_boundaries(md,options_structure,width,i);
    25                 return;
    26         case 'elementnumbering',
    27                 plot_elementnumbering(md,options_structure,width,i);
    28                 return;
    29         case 'highlightelements',
    30                 plot_highlightelements(md,options_structure,width,i);
    31         case 'segmentnumbering',
    32                 plot_segmentnumbering(md,options_structure,width,i);
    33                 return;
    34         case 'importancefactors',
    35                 plot_importancefactors(md,options_structure,width,i);
    36                 return;
    37         case 'elements_type',
    38                 plot_elementstype(md,options_structure,width,i);
    39                 return;
    40         case 'gridnumbering',
    41                 plot_gridnumbering(md,options_structure,width,i);
    42                 return;
    43         case 'highlightgrids',
    44                 plot_highlightgrids(md,options_structure,width,i);
    45                 return;
    46         case 'basal_drag',
    47                 plot_basaldrag(md,options_structure,width,i);
    48                 return;
    49         case 'driving_stress',
    50                 plot_drivingstress(md,options_structure,width,i);
    51                 return;
    52         case 'mesh',
    53                 plot_mesh(md,options_structure,width,i);
    54                 return;
    55         case 'penalties',
    56                 plot_penalties(md,options_structure,width,i);
    57                 return;
    58         case 'quiver',
    59                 plot_quiver(md,options_structure,width,i);
    60                 return;
    61         case 'quiver3',
    62                 plot_quiver3(md,options_structure,width,i);
    63                 return;
    64         case 'quivervel',
    65                 plot_quivervel(md,options_structure,width,i);
    66                 return;
    67         case 'riftvel',
    68                 plot_riftvel(md,options_structure,width,i);
    69                 return;
    70         case 'riftrelvel',
    71                 plot_riftrelvel(md,options_structure,width,i);
    72                 return;
    73         case 'riftpenetration',
    74                 plot_riftpenetration(md,options_structure,width,i);
    75                 return;
    76         case 'sarpwr',
    77                 plot_sarpwr(md,options_structure,width,i)
    78                 return
    79         case {'segmentonneumann_diag','segmentonneumann_prog'}
    80                 plot_segmentonneumann(md,options_structure,width,i,data)
    81                 return
    82         case {'strainrate_tensor','strainrate','strainrate_principal','strainrate_principalaxis1','strainrate_principalaxis2','strainrate_principalaxis3',...
    83                 'stress_tensor','stress','stress_principal','stress_principalaxis1','stress_principalaxis2','stress_principalaxis3',...
    84                 'deviatoricstress_tensor','deviatoricstress','deviatoricstress_principal','deviatoricstress_principalaxis1','deviatoricstress_principalaxis2','deviatoricstress_principalaxis3'},
    85                 plot_tensor(md,options_structure,width,i,data);
    86                 return;
    87         case 'thermaltransient_results',
    88                 plot_thermaltransient_results(md,options_structure,width,i);
    89                 return;
    90         case 'transient_movie',
    91                 plot_transient_movie(md,options_structure,width,i);
    92                 return;
    93         case 'transient_results',
    94                 plot_transient_results(md,options_structure,width,i);
    95                 return;
    96         otherwise,
    97                 if isfield(struct(md),data)
    98                         data=eval(['md.' data ';']);
    99                 else
    100                         error('plot error message: data provided not supported yet. Type plotdoc for help');
    101                 end
    102         end
    103 end
    104 
    105 %Figure out if this is a Section plot
    106 if ~isnan(options_structure.sectionvalue)
    107         plot_section(md,data,options_structure,width,i);
    108         return;
    109 end
    110 
    111 %check on arguments
    112 if (iscell(data) | isempty(data)),
    113         error('plot error message: data provided is empty');
    114 end
    115 
    116 %process data and model
    117 [x y z elements is2d]=processmesh(md,options_structure);
    118 [data isongrid]=processdata(md,data,options_structure);
     7%   See also: PLOTMODEL, PLOT_DEALER
    1198
    1209%edgecolor?
    … …  
    12413        edgecolor='none';
    12514end
    126 
    127 %standard plot:
    128 subplot(width,width,i);
    12915
    13016%element data
    … …  
    14531        if is2d
    14632                A=elements(:,1); B=elements(:,2); C=elements(:,3);
    147                 patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor);
     33                patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', data(:),'FaceColor','interp','EdgeColor',edgecolor);
    14834        else
    14935                if options_structure.layer>=1,
    15036                        A=elements(:,1); B=elements(:,2); C=elements(:,3);
    151                         patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor);
     37                        patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', data(:),'FaceColor','interp','EdgeColor',edgecolor);
    15238                else
    15339                        A=elements(:,1); B=elements(:,2); C=elements(:,3); D=elements(:,4); E=elements(:,5); F=elements(:,6);
    154                         patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor);
    155                         patch( 'Faces', [D E F], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor);
    156                         patch( 'Faces', [A B E D], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor);
    157                         patch( 'Faces', [B E F C ], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor);
    158                         patch( 'Faces', [C A D F ], 'Vertices', [x y z],'FaceVertexCData', data,'FaceColor','interp','EdgeColor',edgecolor);
     40                        patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', data(:),'FaceColor','interp','EdgeColor',edgecolor);
     41                        patch( 'Faces', [D E F], 'Vertices', [x y z],'FaceVertexCData', data(:),'FaceColor','interp','EdgeColor',edgecolor);
     42                        patch( 'Faces', [A B E D], 'Vertices', [x y z],'FaceVertexCData', data(:),'FaceColor','interp','EdgeColor',edgecolor);
     43                        patch( 'Faces', [B E F C ], 'Vertices', [x y z],'FaceVertexCData', data(:),'FaceColor','interp','EdgeColor',edgecolor);
     44                        patch( 'Faces', [C A D F ], 'Vertices', [x y z],'FaceVertexCData', data(:),'FaceColor','interp','EdgeColor',edgecolor);
    15945                end
    16046        end
    16147end
    162 
    163 %apply all options
    164 if isnan(options_structure.shading) & isnan(options_structure.edgecolor) & size(data,1)==md.numberofgrids,
    165         options_structure.shading='interp';
    166 end
    167 
    168 applyoptions(md,data,options_structure);

  • issm/trunk/src/m/classes/public/plot/plotdoc.m

    r1 r27  
    3030disp('                  - ''riftpenetration'': penetration levels for a fault');
    3131disp('                  - ''basal_drag'': plot the basal drag on the bed (in kPa)');
     32disp('                  - ''basal_dragx'' or ''basal_dragy'' : plot a component of the basal drag on the bed (in kPa)');
    3233disp('                  - ''driving_stress'': plot the driving stress (in kPa)');
    3334disp('                  - ''strainrate_tensor'': plot the components of the strainrate tensor (exx,eyy,ezz,exy,exz,eyz) if computed');
    … …  
    5354disp('       ''contourticks'': ''on'' or ''off'' to display the ticks of the contours');
    5455disp('       ''contouronly'': ''on'' or ''off'' to display the contours on a white background');
    55 disp('       ''contourcolors'': ticks and contour color');
     56disp('       ''contourcolor'': ticks and contour color');
     57disp('       ''density'': density of quivers (one arrow every N nodes, N integer)');
     58disp('       ''streamlines'': N (number of stream lines) or {[x1 y1],...} (coordinates of seed points) add streanlines on current figure');
    5659disp('       ''wrapping'': repeat ''n'' times the colormap (''n'' must be an integer)');
    5760disp('       ''edgecolor'': same as standard matlab option EdgeColor (color name: ''black'' or RGB array: [0.5 0.2 0.8])');
    … …  
    6164disp('       ''resolution'': resolution used by section value (array of type [horizontal_resolution vertical_resolution])');
    6265disp('                       horizontal_resolution must be in meter, and vertical_resolution a number of layers');
    63 disp('       ''showsection'': show section used by ''sectionvalue'' (string ''yes'')');
     66disp('       ''showsection'': show section used by ''sectionvalue'' (string ''on'' or a number of labels)');
    6467disp('       ''sectionvalue'': give the value of data on a profile given by an Argus file (string ''Argusfile_name.exp'')');
    6568disp('       ''smooth'': smooth element data (string ''yes'' or integer)');

  • issm/trunk/src/m/classes/public/plot/plotmodel.m

    r1 r27  
    11function plotmodel(md,varargin)
    22%At command prompt, type plotdoc for help on documentation.
     3
    34global ISSM_DIR
    45if isempty(ISSM_DIR),
    … …  
    2829        %recover options  for all subplots.
    2930        options=recover_plot_options(md,varargin,numberofplots);
     31       
     32        %Create figure
     33        figure(figurenumber),clf;
    3034
    31         %Create figure
    32         figure(figurenumber); clf;
    33        
    3435        %Go through all data plottable
    3536        for i=1:numberofplots,
    36                
     37
    3738                %call unit plot
    38                 plot_unit(md,options{i},subplotwidth,i);
    39        
     39                plot_dealer(md,options{i},subplotwidth,i);
     40
    4041        end
    4142else

  • issm/trunk/src/m/classes/public/plot/processdata.m

    r1 r27  
    77%   See also: PLOTMODEL, PROCESSMESH
    88
    9 %some checks
    10 if length(data)~=md.numberofgrids & length(data)~=md.numberofelements & length(data)~=md.numberofgrids*6
     9%check format
     10if (iscell(data) | isempty(data)),
     11        error('plotmodel error message: data provided is empty');
     12end
     13
     14%check length
     15if length(data)~=md.numberofgrids & length(data)~=md.numberofelements & length(data)~=md.numberofgrids*6 & (strcmpi(md.type,'2d') | length(data)~=md.numberofgrids2d)
    1116        error('plotmodel error message: data not supproted yet')
    1217end
    1318
    14 %6*grid data
     19%treat the case length(data)=6*grids
    1520if length(data)==6*md.numberofgrids
    1621        %keep the only norm of data
    … …  
    1823        data2=data(2:6:md.numberofgrids*6);
    1924        data=sqrt(data1.^2+data2.^2);
     25        %---> go to grid data
     26end
     27
     28%treat the case length(data)=grids2d
     29if (strcmpi(md.type,'3d') & length(data)==md.numberofgrids2d),
     30        data=project3d(md,data,'node');
    2031        %---> go to grid data
    2132end
    … …  
    6071        data=project2d(md,data,options_structure.layer); %project onto 2d mesh
    6172end
    62 
    63 

  • issm/trunk/src/m/classes/public/solve.m

    r1 r27  
    1414%some checks on list of arguments
    1515
    16 solutions={'mesh2grid','mesh','thermaltransient','thermalsteady','qmu','diagnostic','diagnostic_horiz','prognostic','transient','paramters','control'};
     16solutions={'mesh2grid','mesh','thermaltransient','thermalsteady','qmu','diagnostic','diagnostic_horiz','prognostic','transient','parameters','control'};
    1717found=0;
    1818for i=1:length(solutions),
    … …  
    2424
    2525if found==0,
    26         error(['solve error messae: solution type ' solutiontype ' not supported yet!']);
     26        error(['solve error message: solution type ' solutiontype ' not supported yet!']);
    2727end
    2828

  • issm/trunk/src/m/classes/public/thicknessevolution.m

    r1 r27  
    44%   This routine compute the new thickness of a model after a time step
    55%   according to the following formula:
    6 %   dh/dt=div(Hu)+Ms-Mb
     6%   dh/dt=div(Hu)
    77%
    88%   Usage:
    … …  
    1313end
    1414
     15%load some variables (it is much faster if the variab;es are loaded from md once for all)
     16numberofelements=md.numberofelements;
     17H=md.thickness;
     18vx=md.vx;
     19vy=md.vy;
     20index=md.elements;
     21x=md.x; y=md.y;
     22
     23%initialization
    1524alpha=zeros(md.numberofelements,3);
    1625beta=zeros(md.numberofelements,3);
     26gradx=zeros(md.numberofgrids,1);
     27grady=zeros(md.numberofgrids,1);
    1728
    18 for n=1:md.numberofelements
    19         X=inv([md.x(md.elements(n,:)) md.y(md.elements(n,:)) ones(3,1)]);
    20         alpha(n,:)=X(1,:);
    21         beta(n,:)=X(2,:);
    22 end
     29%build some usefull variables
     30line=index(:);
     31summation=1/3*ones(3,1);
     32linesize=3*numberofelements;
     33x1=x(index(:,1)); x2=x(index(:,2)); x3=x(index(:,3)); y1=y(index(:,1)); y2=y(index(:,2)); y3=y(index(:,3));
    2334
    24 summation=[1;1;1];
    25 dHu=(md.vx(md.elements).*md.thickness(md.elements).*alpha)*summation;
    26 dHv=(md.vy(md.elements).*md.thickness(md.elements).*beta)*summation;
     35%compute nodal functions coefficients N(x,y)=alpha x + beta y + gamma
     36invdet=1./(x1.*(y2-y3)-x2.*(y1-y3)+x3.*(y1-y2));
     37alpha=[invdet.*(y2-y3) invdet.*(y3-y1) invdet.*(y1-y2)];
     38beta=[invdet.*(x3-x2) invdet.*(x1-x3) invdet.*(x2-x1)];
    2739
    28 dHu=averaging(md,dHu,0);
    29 dHv=averaging(md,dHv,0);
    30 
    31 dhdt=dHu+dHv+md.accumulation-md.melting;
     40%compute dhdt=div(Hu)
     41Hvx=H.*vx; Hvy=H.*vy;
     42dhdt=sum(Hvx(index).*alpha,2)+sum(Hvy(index).*beta,2);
Note: See TracChangeset for help on using the changeset viewer.