Changeset 20548


Ignore:
Timestamp:
04/25/16 10:08:11 (9 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 20546

Location:
issm/trunk
Files:
3 deleted
47 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk

  • issm/trunk/examples/Greenland/runme.m

    r20500 r20548  
    11clear all;
    2 steps=[1 2 3 4 5];
     2steps=[1];
    33
    44%Location of SeaRISE dataset
    … …  
    146146        lat  = ncread(ncbox,'lat');
    147147        lon  = ncread(ncbox,'lon');
    148         smbbox = ncread(ncbox,'SMB');
     148        smbbox = ncread(ncbox,'MassFlux');
    149149        [x1 y1]=ll2xy(lat,lon,+1,45,70);
    150150

  • issm/trunk/examples/ISMIP/CheatyRunme.m

    r19105 r20548  
    33%step 8 is specific to ISMIPF
    44
    5 steps=[1:7];
     5steps=[1];
    66
    77% parameter file to be used, choose between CheatyIsmipA.par or CheatyIsmipF.par

  • issm/trunk/examples/IceBridge/Greenland.par

    r19171 r20548  
    6262
    6363disp('   Interpolating surface mass balance');
    64 md.surfaceforcings.mass_balance=InterpFromGridToMesh(x1,y1,smb,md.mesh.x,md.mesh.y,0);
    65 md.surfaceforcings.mass_balance=md.surfaceforcings.mass_balance*md.materials.rho_water/md.materials.rho_ice;
     64md.smb.mass_balance=InterpFromGridToMesh(x1,y1,smb,md.mesh.x,md.mesh.y,0);
     65md.smb.mass_balance=md.smb.mass_balance*md.materials.rho_water/md.materials.rho_ice;
    6666
    6767disp('   Construct basal friction parameters');

  • issm/trunk/examples/IceBridge/Greenland_cheatsheet.par

    r19171 r20548  
    6666
    6767disp('   Interpolating surface mass balance');
    68 md.surfaceforcings.mass_balance=InterpFromGridToMesh(x1,y1,smb,md.mesh.x,md.mesh.y,0);
    69 md.surfaceforcings.mass_balance=md.surfaceforcings.mass_balance*md.materials.rho_water/md.materials.rho_ice;
     68md.smb.mass_balance=InterpFromGridToMesh(x1,y1,smb,md.mesh.x,md.mesh.y,0);
     69md.smb.mass_balance=md.smb.mass_balance*md.materials.rho_water/md.materials.rho_ice;
    7070
    7171disp('   Construct basal friction parameters');

  • issm/trunk/examples/IceBridge/Greenland_noOIB.par

    r19171 r20548  
    4040
    4141disp('   Interpolating surface mass balance');
    42 md.surfaceforcings.mass_balance=InterpFromGridToMesh(x1,y1,smb,md.mesh.x,md.mesh.y,0);
    43 md.surfaceforcings.mass_balance=md.surfaceforcings.mass_balance*md.materials.rho_water/md.materials.rho_ice;
     42md.smb.mass_balance=InterpFromGridToMesh(x1,y1,smb,md.mesh.x,md.mesh.y,0);
     43md.smb.mass_balance=md.smb.mass_balance*md.materials.rho_water/md.materials.rho_ice;
    4444
    4545disp('   Construct basal friction parameters');

  • issm/trunk/examples/IceBridge/runme.m

    r18267 r20548  
    4646
    4747        md = setmask(md,'','');
    48         md = parameterize(md,'./Greenland_cheatsheet.par');
     48        md = parameterize(md,'./Greenland.par');
    4949        %md = parameterize(md,'./Greenland_noOIB.par');
    5050        md = setflowequation(md,'SSA','all');
    … …  
    108108        smb = InterpFromGridToMesh(x1,y1,smb',md.mesh.x,md.mesh.y,0)*1000/md.materials.rho_ice;
    109109        smb = [smb smb smb-1.0];
    110         md.surfaceforcings.mass_balance = [smb;1 10 20];
     110        md.smb.mass_balance = [smb;1 10 20];
    111111
    112112        %Set transient options, run for 20 years, saving every year
    … …  
    117117        %Additional options
    118118        md.inversion.iscontrol=0;
    119         md.transient.requested_outputs={'IceVolume','TotalSmb','SurfaceforcingsMassBalance'};
     119        md.transient.requested_outputs={'IceVolume','TotalSmb','SmbMassBalance'};
    120120        md.verbose=verbose('solution',true,'module',true,'convergence',true);
    121121
    … …  
    134134        plotmodel(md,'data',md.results.TransientSolution(end).Vel,'caxis',[1e-1 6000],...
    135135                'log', 10, 'title', 'Velocity (m/y)','gridded',1, ...
    136                 'data', md.results.TransientSolution(end).SurfaceforcingsMassBalance, ...
     136                'data', md.results.TransientSolution(end).SmbMassBalance, ...
    137137                'title', 'Surface mass balance (m/y)','gridded',1, ...
    138138                'data',md.results.TransientSolution(end).Thickness,...
    … …  
    146146        %Plot surface mass balance
    147147        surfmb=[]; for i=1:100; surfmb=[surfmb ...
    148                 md.results.TransientSolution(i).SurfaceforcingsMassBalance]; end
     148                md.results.TransientSolution(i).SmbMassBalance]; end
    149149        subplot(3,1,1); plot([0.2:0.2:20],mean(surfmb)); title('Mean Surface mass balance');
    150150

  • issm/trunk/examples/IceflowModels/EISMINT.par

    r18198 r20548  
    2929sb=10^-2/1000.; %m/yr/m
    3030rel=450.*1000.; %m
    31 md.surfaceforcings.mass_balance=min(smb_max,sb*(rel-radius));
     31md.smb.mass_balance=min(smb_max,sb*(rel-radius));
    3232
    3333disp('      creating velocities');

  • issm/trunk/examples/Jakobshavn/Jks.par

    r19171 r20548  
    3838
    3939disp('   Interpolating surface mass balance');
    40 md.surfaceforcings.mass_balance=InterpFromGridToMesh(x1,y1,smb,md.mesh.x,md.mesh.y,0);
    41 md.surfaceforcings.mass_balance=md.surfaceforcings.mass_balance*md.materials.rho_water/md.materials.rho_ice;
     40md.smb.mass_balance=InterpFromGridToMesh(x1,y1,smb,md.mesh.x,md.mesh.y,0);
     41md.smb.mass_balance=md.smb.mass_balance*md.materials.rho_water/md.materials.rho_ice;
    4242
    4343disp('   Construct basal friction parameters');

  • issm/trunk/examples/Jakobshavn/runme.m

    r18301 r20548  
    1 steps=[4];
     1steps=[1];
    22
    33if any(steps==1)

  • issm/trunk/examples/UncertaintyQuantification/runme.m

    r18269 r20548  
    11%PIG Uncertainty Quantification Application
    2 steps=[1:3];
     2steps=[1];
    33
    44if any(steps==1)
    … …  
    140140
    141141        md.stressbalance.restol=10^-5; %tighten tolerances for UQ analyses
     142
     143        %Turn off verbosity
     144        md.verbose=verbose(0);
    142145       
    143146        %solve
    144147        md.qmu.isdakota=1; md.inversion.iscontrol=0;
    145148        md.cluster=generic('name',oshostname,'np',4);
     149
     150        %Dakota runs in parallel with a master/slave configuration.
     151        % At least 2 cpu's are needed to run the UQ
     152        md.qmu.params.evaluation_scheduling='master';
     153        md.qmu.params.processors_per_evaluation=md.cluster.np-1;
     154
    146155        md=solve(md,StressbalanceSolutionEnum,'overwrite','y');
    147156
    … …  
    217226        md.inversion.iscontrol=0;
    218227        md.cluster=generic('name',oshostname,'np',2);
     228
     229        %Dakota runs in parallel with a master/slave configuration.
     230        % At least 2 cpu's are needed to run the UQ
     231        md.qmu.params.evaluation_scheduling='master';
     232        md.qmu.params.processors_per_evaluation=md.cluster.np-1
     233
    219234        md.verbose=verbose('qmu',true);
    220235        md=solve(md,StressbalanceSolutionEnum,'overwrite','y');

  • issm/trunk/jenkins/jenkins.sh

    r20507 r20548  
    77rm -rf $ISSM_DIR/nightlylog
    88mkdir  $ISSM_DIR/nightlylog
     9
     10#Server URI
     11SERVER='https://ross.ics.uci.edu:8080'
    912
    1013#Get configuration
    … …  
    4447        #Get changes from jenkins itself (svn requires credentials)
    4548        rm -rf changes
    46         wget http://ross.ics.uci.edu:8080/job/$JOB_NAME/$BUILD_NUMBER/changes > /dev/null 2>&1
     49        wget $SERVER/job/$JOB_NAME/$BUILD_NUMBER/changes > /dev/null 2>&1
    4750
    4851        #Process html page and get the list of files that has changed (tricky...)
    … …  
    6063        if [ ! -z "$(cat $ISSM_DIR/TEMP | grep externalpackages)" ] ; then
    6164                echo "  -- checking for changed externalpackages... yes";
    62                 ISSM_EXTERNALPACKAGES="install"
     65                ISSM_EXTERNALPACKAGES="yes"
    6366        else
    6467                echo "  -- checking for changed externalpackages... no";
    65                 ISSM_EXTERNALPACKAGES="none"
     68                ISSM_EXTERNALPACKAGES="no"
    6669        fi
    6770
    … …  
    8992        echo "  -- checking for reconfiguration... yes";
    9093        echo "  -- checking for recompilation... yes";
    91         ISSM_EXTERNALPACKAGES="install"
     94        ISSM_EXTERNALPACKAGES="yes"
    9295        ISSM_RECONFIGURE="yes"
    9396        ISSM_COMPILATION="yes"
    … …  
    97100echo $SVN_REVISION_1 > $ISSM_DIR/svn_revision_old
    98101#}}}
    99 #install/none external packages    (ISSM_EXTERNALPACKAGES){{{
     102#Install external packages    (ISSM_EXTERNALPACKAGES){{{
    100103
    101104#Jenkins xml files for individual packages
    … …  
    115118        cd $ISSM_DIR/externalpackages/$PACKAGENAME
    116119
    117         #install if requested or if install does not exist
    118         if [ "$ISSM_EXTERNALPACKAGES" == "install" ] || [ ! -d install ] || [ ! "$(ls -A install)" ];
    119         then
     120        #install if requested or if previous install has not been successful
     121        if [ "$ISSM_EXTERNALPACKAGES" == "yes" ] || [ ! -f SUCCESS ]; then
    120122                echo "======================================================";
    121123                echo "       Installing $PACKAGENAME                        ";
    … …  
    126128                        echo "    ERROR: installation of $PACKAGENAME failed        ";
    127129                        echo "======================================================";
    128                         #erase install directory, so that next time, we still try and compile this!
    129                         rm -rf install
    130130                        echo "<testcase classname=\"externalpackages\" name=\"$PACKAGENAME\">" >> $EXTERNAL_TEST_FILE
    131131                        echo '<failure message="failure">External packages did not install right. Check it.</failure>' >> $EXTERNAL_TEST_FILE
    132132                        echo '</testcase>' >> $EXTERNAL_TEST_FILE
     133                        exit 1
    133134                else
    134135                        echo "<testcase classname=\"externalpackages\" name=\"$PACKAGENAME\"/>" >> $EXTERNAL_TEST_FILE
     136                        touch SUCCESS
    135137                fi
    136138                source $ISSM_DIR/etc/environment.sh

  • issm/trunk/src

  • issm/trunk/src/c/Makefile.am

    r20500 r20548  
    534534                        ./shared/TriMesh/SplitMeshForRifts.cpp\
    535535                        ./shared/TriMesh/TriMeshUtils.cpp\
    536                         ./modules/TriaSearchx/TriaSearchx.cpp\
    537536                        ./modules/TriMeshx/TriMeshx.cpp\
    538537                        ./modules/TriMeshProcessRiftsx/TriMeshProcessRiftsx.cpp\
    … …  
    594593if !STANDALONE_LIBRARIES
    595594libISSMCore_la_LIBADD = $(PETSCLIB) $(TAOLIB) $(M1QN3LIB) $(PLAPACKLIB) $(MUMPSLIB) $(SUPERLULIB) $(SPOOLESLIB) $(SCALAPACKLIB) $(BLACSLIB) $(HYPRELIB) $(SPAILIB) $(PROMETHEUSLIB) $(PASTIXLIB) $(MLLIB) $(DAKOTALIB) $(METISLIB) $(CHACOLIB) $(SCOTCHLIB) $(BLASLAPACKLIB) $(MKLLIB) $(MPILIB) $(MATHLIB) $(GRAPHICSLIB) $(MULTITHREADINGLIB) $(OSLIBS) $(GSLLIB)   $(ADOLCLIB) $(AMPILIB) $(METEOIOLIB) $(SNOWPACKLIB)
     595if FORTRAN
     596libISSMCore_la_LIBADD += $(FLIBS) $(FORTRANLIB)
     597endif
    596598endif
    597599endif
    … …  
    645647libISSMOverload_la_CFLAGS  = -D_C_ -fPIC $(COPTFLAGS) $(CFLAGS)
    646648#}}}
    647 
    648649#Executable {{{
    649650if ANDROID

  • issm/trunk/src/c/analyses/LevelsetAnalysis.cpp

    r20500 r20548  
    7777/*Finite element Analysis*/
    7878void           LevelsetAnalysis::Core(FemModel* femmodel){/*{{{*/
    79 
    80         #if !defined(_DEVELOPMENT_)
    81         _error_("Not implemented yet");
    82         #endif
    8379
    8480        /*parameters: */

  • issm/trunk/src/c/classes/Elements/PentaRef.cpp

    r19105 r20548  
    443443
    444444        /*Get nodal functions derivatives in reference triangle*/
    445         IssmDouble* dbasis_ref=xNew<IssmDouble>(3*numnodes);
     445        IssmDouble dbasis_ref[3*NUMNODESMAX];
    446446        GetNodalFunctionsDerivativesReference(dbasis_ref,gauss,finiteelement);
    447447
    … …  
    461461                dbasis[numnodes*2+i]=Jinv[2][0]*dbasis_ref[0*numnodes+i]+Jinv[2][1]*dbasis_ref[1*numnodes+i]+Jinv[2][2]*dbasis_ref[2*numnodes+i];
    462462        }
    463 
    464         /*Clean up*/
    465         xDelete<IssmDouble>(dbasis_ref);
    466463}
    467464/*}}}*/

  • issm/trunk/src/c/classes/Elements/SegRef.cpp

    r19105 r20548  
    146146
    147147        /*Get nodal functions derivatives in reference triangle*/
    148         IssmDouble* dbasis_ref=xNew<IssmDouble>(numnodes);
     148        IssmDouble dbasis_ref[1*NUMNODESMAX];
    149149        GetNodalFunctionsDerivativesReference(dbasis_ref,gauss,finiteelement);
    150150
    … …  
    158158                dbasis[i] = Jinv*dbasis_ref[i];
    159159        }
    160 
    161         /*Clean up*/
    162         xDelete<IssmDouble>(dbasis_ref);
    163 
    164160}
    165161/*}}}*/

  • issm/trunk/src/c/classes/Elements/TetraRef.cpp

    r19105 r20548  
    225225
    226226        /*Get nodal functions derivatives in reference triangle*/
    227         IssmDouble* dbasis_ref=xNew<IssmDouble>(3*numnodes);
     227        IssmDouble dbasis_ref[3*NUMNODESMAX];
    228228        GetNodalFunctionsDerivativesReference(dbasis_ref,gauss,finiteelement);
    229229
    … …  
    243243                dbasis[numnodes*2+i]=Jinv[2][0]*dbasis_ref[0*numnodes+i]+Jinv[2][1]*dbasis_ref[1*numnodes+i]+Jinv[2][2]*dbasis_ref[2*numnodes+i];
    244244        }
    245 
    246         /*Clean up*/
    247         xDelete<IssmDouble>(dbasis_ref);
    248245}
    249246/*}}}*/

  • issm/trunk/src/c/classes/Elements/TriaRef.cpp

    r19105 r20548  
    197197
    198198        /*Get nodal functions derivatives in reference triangle*/
    199         IssmDouble* dbasis_ref=xNew<IssmDouble>(2*numnodes);
     199        IssmDouble dbasis_ref[2*NUMNODESMAX];
    200200        GetNodalFunctionsDerivativesReference(dbasis_ref,gauss,finiteelement);
    201201
    … …  
    211211                dbasis[numnodes*1+i] = Jinv[1][0]*dbasis_ref[0*numnodes+i]+Jinv[1][1]*dbasis_ref[1*numnodes+i];
    212212        }
    213 
    214         /*Clean up*/
    215         xDelete<IssmDouble>(dbasis_ref);
    216213
    217214}

  • issm/trunk/src/c/modules/modules.h

    r20500 r20548  
    8989#include "./SpcNodesx/SpcNodesx.h"
    9090#include "./SurfaceAreax/SurfaceAreax.h"
    91 #include "./TriaSearchx/TriaSearchx.h"
    9291#include "./TriMeshx/TriMeshx.h"
    9392#include "./TriMeshProcessRiftsx/TriMeshProcessRiftsx.h"

  • issm/trunk/src/m/boundaryconditions/SetMarineIceSheetBC.m

    r20500 r20548  
    4040pos=find(md.mesh.vertexonboundary & ~vertexonicefront);
    4141if isempty(pos),
    42         warning('SetMarineIceSheetBC warning: ice front all around the glacier, no dirichlet found. Dirichlet must be added manually')
     42        disp('SetMarineIceSheetBC warning: ice front all around the glacier, no dirichlet applied')
    4343end
    4444md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1);

  • issm/trunk/src/m/classes/calving.m

    r20500 r20548  
    3939                        if (solution~=TransientSolutionEnum() | md.transient.ismovingfront==0), return; end
    4040
    41                         md = checkfield(md,'fieldname','calving.calvingrate(1:md.mesh.numberofvertices,:)','>=',0,'timeseries',1,'NaN',1,'Inf',1);
    42                         md = checkfield(md,'fieldname','calving.meltingrate(1:md.mesh.numberofvertices,:)','>=',0,'timeseries',1,'NaN',1,'Inf',1);
     41                        md = checkfield(md,'fieldname','calving.calvingrate','>=',0,'timeseries',1,'NaN',1,'Inf',1);
     42                        md = checkfield(md,'fieldname','calving.meltingrate','>=',0,'timeseries',1,'NaN',1,'Inf',1);
    4343                end % }}}
    4444                function disp(self) % {{{

  • issm/trunk/src/m/classes/clusters/generic_static.m

    r20500 r20548  
    5151                function BuildQueueScript(cluster,dirname,modelname,solution,io_gather,isvalgrind,isgprof,isdakota) % {{{
    5252
     53                        %Check that issm.exe exists in the right path
     54                        if ~exist([cluster.codepath '/issm.exe'],'file'),
     55                                error(['File ' cluster.codepath '/issm.exe does not exist']);
     56                        end
     57
     58                        %Now process codepath and replace empty spaces with \ to avoid errors in queuing script
     59                        codepath2=strrep(cluster.codepath,' ','\ ');
     60
     61                        %write queuing script
     62                        %what is the executable being called?
     63                        executable='issm.exe';
     64                        if isdakota,
     65                                version=IssmConfig('_DAKOTA_VERSION_'); version=str2num(version(1:3));
     66                                if (version>=6),
     67                                        executable='issm_dakota.exe';
     68                                end
     69                        end
     70
    5371                        %write queuing script
    5472                        fid=fopen([modelname '.queue'],'w');
    5573                        fprintf(fid,'#!%s\n',cluster.shell);
    56                         fprintf(fid,[cluster.codepath '/mpiexec -np %i %s/issm.exe %s %s %s \n'],cluster.np,cluster.codepath,EnumToString(solution),'./',modelname);
     74                        fprintf(fid,['%s/mpiexec -np %i %s/%s %s %s %s \n'],codepath2,cluster.np,codepath2,executable,EnumToString(solution),'./',modelname);
    5775                        fclose(fid);
    5876

  • issm/trunk/src/m/classes/clusters/lonestar.m

    r20500 r20548  
    123123                                 fprintf(fid,['module load ' cluster.modules{i} '\n']);
    124124                         end
     125
     126                         %fprintf(fid,'export KMP_AFFINITY="granularity=fine,compact,verbose" \n\n');
     127
    125128                         if length(find(cluster.email=='@'))>0
    126                                  fprintf(fid,'#SBATCH --mail-user=%s \n',cluster.email);
    127                                  fprintf(fid,'#SBATCH --mail-type=end \n\n');
     129                                 %fprintf(fid,'#SBATCH --mail-user=%s \n',cluster.email);
     130                                 %fprintf(fid,'#SBATCH --mail-type=end \n\n');
     131
     132                                 fprintf(fid,'ssh login1 "mail -s ''SLURM Jobid=${SLURM_JOBID} Name=${SLURM_JOB_NAME} Began on Lonestar 5.'' %s <<< ''Job Started'' " \n\n',cluster.email);
    128133                         end
     134
    129135                         fprintf(fid,'export PATH="$PATH:."\n\n');
    130136                         fprintf(fid,'export ISSM_DIR="%s/../"\n',cluster.codepath); %FIXME
    … …  
    134140                         if ~io_gather, %concatenate the output files:
    135141                                 fprintf(fid,'cat %s.outbin.* > %s.outbin',modelname,modelname);
     142                         end
     143
     144                         if length(find(cluster.email=='@'))>0
     145                                 fprintf(fid,'\n');
     146                                 fprintf(fid,'ssh login1 "mail -s ''SLURM Jobid=${SLURM_JOBID} Name=${SLURM_JOB_NAME} Ended on Lonestar 5.'' %s <<< ''Job Ended'' " \n\n',cluster.email);
    136147                         end
    137148                         fclose(fid);

  • issm/trunk/src/m/classes/clusters/pfe.m

    r20500 r20548  
    9292                                         end
    9393                                 end
    94                          
     94                         elseif strcmpi(cluster.processor,'bro'),
     95                                 if cluster.hyperthreading,
     96                                         if ((cluster.cpuspernode>56 ) | (cluster.cpuspernode<1)),
     97                                                 md = checkmessage(md,'cpuspernode should be between 1 and 56 for ''bro'' processors in hyperthreading mode');
     98                                         end
     99                                 else
     100                                         if ((cluster.cpuspernode>28 ) | (cluster.cpuspernode<1)),
     101                                                 md = checkmessage(md,'cpuspernode should be between 1 and 28 for ''bro'' processors');
     102                                         end
     103                                 end
    95104                         elseif strcmpi(cluster.processor,'has'),
    96105                                 if cluster.hyperthreading,

  • issm/trunk/src/m/classes/independent.py

    r20500 r20548  
    22from pairoptions import pairoptions
    33from fielddisplay import fielddisplay
     4from checkfield import checkfield
    45from MatlabFuncs import *
    56from EnumDefinitions import *
    … …  
    5253                                raise TypeError("independent checkconsistency error: nods should be set to the size of the independent variable")
    5354
    54                 if self.fov_forward_indices:
     55                if len(self.fov_forward_indices) > 0:
    5556                        if not strcmpi(driver,'fov_forward'):
    5657                                raise TypeError("cannot declare an independent with fov_forward_indices when the driver is not fov_forward!")

  • issm/trunk/src/m/classes/miscellaneous.m

    r20500 r20548  
    3939                end % }}}
    4040                function savemodeljs(self,fid,modelname) % {{{
    41                
     41
    4242                        writejsstring(fid,[modelname '.miscellaneous.notes'],self.notes);
    4343                        writejsstring(fid,[modelname '.miscellaneous.name'],self.name);
    4444                        if strcmpi(class(self.dummy),'double'),
    45                                 writejs1Darray(fid,[modelname '.miscellaneous.dummy'],self.dummy);
     45                                if size(self.dummy,2)==1,
     46                                        writejs1Darray(fid,[modelname '.miscellaneous.dummy'],self.dummy);
     47                                else
     48                                        writejs2Darray(fid,[modelname '.miscellaneous.dummy'],self.dummy);
     49                                end
    4650                        end
    4751

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

    r20500 r20548  
    208208                                md.smb.smbref=project2d(md,md.smb.smbref,md.mesh.numberoflayers);
    209209                        end;
    210                         if ~isnan(md.balancethickness.thickening_rate), md.balancethickness.thickening_rate=project2d(md,md.balancethickness.thickening_rate,md.mesh.numberoflayers); end;
    211210
    212211                        %results

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

    r20500 r20548  
    687687                        md.smb.mass_balance=project2d(md,md.smb.mass_balance,md.mesh.numberoflayers)
    688688               
    689                 if not numpy.isnan(md.balancethickness.thickening_rate).all(): md.balancethickness.thickening_rate=project2d(md,md.balancethickness.thickening_rate,md.mesh.numberoflayers)
    690 
    691689                #results
    692690                if not numpy.isnan(md.initialization.vx).all(): md.initialization.vx=DepthAverage(md,md.initialization.vx)

  • issm/trunk/src/m/coordsystems/gdaltransform.m

    r20500 r20548  
    66%
    77%   Example:
    8 %      [x,y] = gdaltransform(md.mesh.lat,md.mesh.long,'EPSG:3184','EPSG:3411');
     8%      [x,y] = gdaltransform(md.mesh.long,md.mesh.lat,'EPSG:4326','EPSG:3031')
    99%
    1010%   For reference:
    11 %       EPSG: 4326 (lat,long)
    12 %       EPSG: 3411  (greenland, +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs)
    13 %       EPSG: 3031 (antarctica, +proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs)
     11%      EPSG: 4326 (lat,long)
     12%      EPSG: 3411 (Greenland,  UPS 45W, 70N)
     13%      EPSG: 3031 (Antarctica, UPS 0E,  71S)
     14%
     15%      ll2xy default projection Antarctica:
     16%        +proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.448564109 +units=m +no_defs
     17%      ll2xy default projection Greenland:
     18%        +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.448564109 +units=m +no_defs
     19%      Bamber's Greeland projection
     20%        +proj=stere +lat_0=90 +lat_ts=71 +lon_0=-39 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.448564109 +units=m +no_defs
     21%
     22%      To get proj.4 string from EPSG, use gdalsrsinfo. Example:    gdalsrsinfo "EPSG:4326" | grep "PROJ.4"
    1423
    1524        %give ourselves a unique temporary directory:
    … …  
    1726
    1827        fid=fopen([temproot '/.rand1234.txt'],'w');
    19         for i=1:length(x),
    20                 fprintf(fid,'%g %g\n',x(i),y(i));
    21         end
     28        fprintf(fid,'%8g %8g\n',[x(:) y(:)]');
    2229        fclose(fid);
    2330
    2431        [s,r]=system(['gdaltransform -s_srs ',proj_in,' -t_srs ',proj_out,'  < ' temproot '/.rand1234.txt > ' temproot '/.rand1235.txt']);
    25         A=textread([temproot '/.rand1235.txt']);
    26         xout=A(:,1);
    27         yout=A(:,2);
     32        if s~=0 | ~isempty(deblank(r)),
     33                error(r);
     34        end
     35        A=load([temproot '/.rand1235.txt']);
     36        xout=A(:,1); xout=reshape(xout,size(x));
     37        yout=A(:,2); yout=reshape(yout,size(y));

  • issm/trunk/src/m/exp/flowlines.m

    r15396 r20548  
    33%
    44%   Usage:
    5 %      flowpath=flowlines(index,x,y,u,v,x0,y0)
     5%      flowpath=flowlines(index,x,y,u,v,x0,y0,varargin)
    66%
    77%   the velocity field is given by the couple (u,v) and the coordinates
    … …  
    1010%
    1111%   Example:
    12 %      flowpath=flowlines(md.mesh.elements,md.mesh.x,md.mesh.y,md.vx,md.initialization.vy,x0,y0)
    13 
    14 %check input size
    15 if nargin>8 | nargin<7,
    16         help flowlines
    17         error('flowlines error message: bad usage');
    18 end
     12%      flowpath=flowlines(md.mesh.elements,md.mesh.x,md.mesh.y,md.initialization.vx,md.initialization.vy,x0,y0)
     13%
     14%   Options:
     15%      - 'maxiter':   how many steps upstream and downstream of the seed points (default: 200)
     16%      - 'precision': division of each segment (higer precision increases number of segments, default: 1)
     17%      - 'downstream':flow line upstream of the seed points (default: 1)
     18%      - 'upstream':  flow line upstream of the seed points (default: 1)
    1919
    2020%check input
    … …  
    2929end
    3030
    31 %get maxiter and precision
    32 if nargin==8
    33         maxiter=varargin{1};
    34 else
    35         maxiter=200; %maximum number of iterations
    36 end
    37 precision=1; %division of each segment (higer precision increases number of segments)
     31%process options
     32options    = pairoptions(varargin{:});
     33maxiter    = getfieldvalue(options,'maxiter',200);
     34precision  = getfieldvalue(options,'precision',1);
     35downstream = getfieldvalue(options,'downstream',1);
     36upstream   = getfieldvalue(options,'upstream',1);
    3837
    39 %check seed points
    40 tria=TriaSearch(index,x,y,x0,y0);
    41 %tria=tsearch(x,y,index,x0,y0);
     38%Create triangulation once for all and check seed points
     39trep = triangulation(index,x,y);
     40tria = pointLocation(trep,[x0 y0]);
    4241pos=find(isnan(tria));
    4342x0(pos)=[];
    … …  
    6362v=v(index)*[1;1;1]/3;
    6463
    65 %initialization:
    66 counter=1;
     64if downstream,
     65        %initialization:
     66        counter=1;
    6767
    68 while any(~done)
     68        while any(~done)
    6969
    70         %find current triangle
    71         queue=find(~done);
    72         tria=TriaSearch(index,x,y,X(queue),Y(queue));
    73         %tria=tsearch(x,y,index,X(queue),Y(queue));
     70                %find current triangle
     71                queue=find(~done);
     72                tria = pointLocation(trep,[X(queue),Y(queue)]);
    7473
    75         %check that the point is actually inside a triangle of the mesh
    76         listnan=find(isnan(tria));
    77         for i=1:length(listnan)
    78                 %remove the last point
    79                 flowpath(queue(listnan(i))).x(end)=[];
    80                 flowpath(queue(listnan(i))).y(end)=[];
    81                 done(queue(listnan(i)))=1;
    82         end
    83         tria(listnan)=[];
    84         queue(listnan)=[];
     74                %check that the point is actually inside a triangle of the mesh
     75                listnan=find(isnan(tria));
     76                for i=1:length(listnan)
     77                        %remove the last point
     78                        flowpath(queue(listnan(i))).x(end)=[];
     79                        flowpath(queue(listnan(i))).y(end)=[];
     80                        done(queue(listnan(i)))=1;
     81                end
     82                tria(listnan)=[];
     83                queue(listnan)=[];
    8584
    86         if isempty(tria),
    87                 break;
    88         end
     85                if isempty(tria),
     86                        break;
     87                end
    8988
    90         %velocity of the current triangle and norm it
    91         ut=u(tria); vt=v(tria); normv=sqrt(ut.^2+vt.^2);
    92         ut=ut./normv;vt=vt./normv;
     89                %velocity of the current triangle and norm it
     90                ut=u(tria); vt=v(tria); normv=max(eps,sqrt(ut.^2+vt.^2));
     91                ut=ut./normv;vt=vt./normv;
    9392
    94         %check counter
    95         if counter>maxiter
    96                 disp(['Maximum number of iterations (' num2str(maxiter) ') reached while going forward'])
    97                 break
    98         end
    99         counter=counter+1;
     93                %check counter
     94                if counter>maxiter
     95                        disp(['Maximum number of iterations (' num2str(maxiter) ') reached while going forward'])
     96                        break
     97                end
     98                counter=counter+1;
    10099
    101         %remove stagnant point
    102         done(queue(find(ut==0 & vt==0)))=1;
     100                %remove stagnant point
     101                done(queue(find(ut==0 & vt==0)))=1;
    103102
    104         %build next point
    105         for i=1:length(queue)
    106                 X(queue(i))=flowpath(queue(i)).x(end)+ut(i)*length_tria(tria(i))/precision;
    107                 Y(queue(i))=flowpath(queue(i)).y(end)+vt(i)*length_tria(tria(i))/precision;
    108                 flowpath(queue(i)).x=[flowpath(queue(i)).x;flowpath(queue(i)).x(end)+ut(i)*length_tria(tria(i))/precision];
    109                 flowpath(queue(i)).y=[flowpath(queue(i)).y;flowpath(queue(i)).y(end)+vt(i)*length_tria(tria(i))/precision];
     103                %build next point
     104                for i=1:length(queue)
     105                        X(queue(i))=flowpath(queue(i)).x(end)+ut(i)*length_tria(tria(i))/precision;
     106                        Y(queue(i))=flowpath(queue(i)).y(end)+vt(i)*length_tria(tria(i))/precision;
     107                        flowpath(queue(i)).x=[flowpath(queue(i)).x;flowpath(queue(i)).x(end)+ut(i)*length_tria(tria(i))/precision];
     108                        flowpath(queue(i)).y=[flowpath(queue(i)).y;flowpath(queue(i)).y(end)+vt(i)*length_tria(tria(i))/precision];
     109                end
    110110        end
    111111end
    112112
    113113%same process but reverse (vel=-vel) to have a vcomplete flow line
    114 queue=[];
    115 counter=1;
    116 X=x0; Y=y0;
    117 done=zeros(N,1);
     114if upstream,
     115        queue=[];
     116        counter=1;
     117        X=x0; Y=y0;
     118        done=zeros(N,1);
    118119
    119 while any(~done)
     120        while any(~done)
    120121
    121         %find current triangle
    122         queue=find(~done);
    123         tria=TriaSearch(index,x,y,X(queue),Y(queue));
    124         %tria=tsearch(x,y,index,X(queue),Y(queue));
     122                %find current triangle
     123                queue=find(~done);
     124                tria = pointLocation(trep,[X(queue),Y(queue)]);
    125125
    126         %check that the point is actually inside a triangle of the mesh
    127         listnan=find(isnan(tria));
    128         for i=1:length(listnan)
    129                 %remove the last point
    130                 flowpath(queue(listnan(i))).x(1)=[];
    131                 flowpath(queue(listnan(i))).y(1)=[];
    132                 done(queue(listnan(i)))=1;
    133         end
    134         tria(listnan)=[];
    135         queue(listnan)=[];
     126                %check that the point is actually inside a triangle of the mesh
     127                listnan=find(isnan(tria));
     128                for i=1:length(listnan)
     129                        %remove the last point
     130                        flowpath(queue(listnan(i))).x(1)=[];
     131                        flowpath(queue(listnan(i))).y(1)=[];
     132                        done(queue(listnan(i)))=1;
     133                end
     134                tria(listnan)=[];
     135                queue(listnan)=[];
    136136
    137         if isempty(tria),
    138                 break;
    139         end
     137                if isempty(tria),
     138                        break;
     139                end
    140140
    141         %velocity of the current triangle and norm it
    142         ut=-u(tria); vt=-v(tria); normv=sqrt(ut.^2+vt.^2);
    143         ut=ut./normv;vt=vt./normv;
     141                %velocity of the current triangle and norm it
     142                ut=-u(tria); vt=-v(tria); normv=max(eps,sqrt(ut.^2+vt.^2));
     143                ut=ut./normv;vt=vt./normv;
    144144
    145         %check counter
    146         if counter>maxiter
    147                 disp(['Maximum number of iterations (' num2str(maxiter) ') reached while going backward'])
    148                 break
    149         end
    150         counter=counter+1;
     145                %check counter
     146                if counter>maxiter
     147                        disp(['Maximum number of iterations (' num2str(maxiter) ') reached while going backward'])
     148                        break
     149                end
     150                counter=counter+1;
    151151
    152         %remove stagnant point
    153         done(queue(find(ut==0 & vt==0)))=1;
     152                %remove stagnant point
     153                done(queue(find(ut==0 & vt==0)))=1;
    154154
    155         %build next point
    156         for i=1:length(queue)
    157                 X(queue(i))=flowpath(queue(i)).x(1)+ut(i)*length_tria(tria(i))/precision;
    158                 Y(queue(i))=flowpath(queue(i)).y(1)+vt(i)*length_tria(tria(i))/precision;
    159                 flowpath(queue(i)).x=[flowpath(queue(i)).x(1)+ut(i)*length_tria(tria(i))/precision; flowpath(queue(i)).x];
    160                 flowpath(queue(i)).y=[flowpath(queue(i)).y(1)+vt(i)*length_tria(tria(i))/precision; flowpath(queue(i)).y];
     155                %build next point
     156                for i=1:length(queue)
     157                        X(queue(i))=flowpath(queue(i)).x(1)+ut(i)*length_tria(tria(i))/precision;
     158                        Y(queue(i))=flowpath(queue(i)).y(1)+vt(i)*length_tria(tria(i))/precision;
     159                        flowpath(queue(i)).x=[flowpath(queue(i)).x(1)+ut(i)*length_tria(tria(i))/precision; flowpath(queue(i)).x];
     160                        flowpath(queue(i)).y=[flowpath(queue(i)).y(1)+vt(i)*length_tria(tria(i))/precision; flowpath(queue(i)).y];
     161                end
    161162        end
    162163end
    163164
    164 %EXP compatibility
     165%EXP compatibility (add name)
    165166for i=1:length(queue)
    166167        flowpath(queue(i)).name=['flowline' num2str(i)];

  • issm/trunk/src/m/mesh/patchglobe.m

    r20500 r20548  
    66        %recover basic options:
    77        bandwidth=getfieldvalue(options,'bandwidth',100000);
     8
     9        %some checks on the mesh:
     10        if (isempty(mh.x) | isempty(mh.y) | isempty(mh.z) | isempty(mh.lat) | isempty(mh.long) | isempty(mh.r) )
     11                error('patchglobe error message: 3D planet mesh has one of the following empty: ''x'',''y'',''z'',''lat'',''long'' or ''r''');
     12        end
     13        if (isempty(mh2d.x) | isempty(mh2d.y) | isempty(mh2d.lat) | isempty(mh2d.long))
     14                error('patchglobe error message: 3D planet mesh has one of the following empty: ''x'',''y'',''lat'' or ''long''');
     15        end
    816
    917        %give ourselves a unique temporary directory:

  • issm/trunk/src/m/miscellaneous/transientrestart.m

    r20500 r20548  
    1818
    1919%Change time
    20 md.timestepping.start_time = results.time/md.constants.yts;
     20md.timestepping.start_time = results.time;
    2121
    2222%Change initialization fields

  • issm/trunk/src/m/plot/applyoptions.js

    r20500 r20548  
    1111                if (options.getfieldvalue('colorbar')==1) {
    1212                        //Handle movie data {{{
    13                         if (typeof data == 'string') {
     13                                if (typeof data == 'string') {
    1414                                //throw Error('plot error message: data of type string');
    1515                                return;
    … …  
    2525                                        })
    2626                                });
    27                                 data = data[node['movieFrame']].slice(0,-1);
     27                                data = data[0].slice(0,-1);
    2828                        } //}}}
    2929                        //Variable options initialization {{{
     30                       
    3031                        var caxis = options.getfieldvalue('caxis',[ArrayMin(data),ArrayMax(data)]);
    3132                        var canvassize = options.getfieldvalue('canvassize',480);
    … …  
    4748                        var clabelitem;
    4849                        var precision=options.getfieldvalue('colorbarprecision',3);
    49                         if (options.exist('log')) {
     50                       
     51                        if (options.getfieldvalue('log',10)!='off') {
    5052                                for (var i = cdivisions; i >= 0; i--) {
    51                                         if (caxisdelta*i/cdivisions==min) {
    52                                                 labels[i] = (caxis[0]).toPrecision(precision);
    53                                         }
    54                                         else {
    55                                                 labels[i] = (Math.exp(Math.log(caxisdelta)*i/cdivisions)+caxis[0]).toPrecision(precision);
    56                                         }       
     53                                        var scale = Math.log(caxis[1])-Math.log(caxis[0]);
     54                                        labels[i] = (Math.exp(Math.log(caxis[0])+scale*(cdivisions-i)/cdivisions)).toPrecision(precision);
    5755                                }
    5856                        } else {
    … …  
    228226        var cmap = options.getfieldvalue('cmap','jet');
    229227        var colorbar = colorbars[cmap];
    230         for (var i = 0; i < colorbar.length; i++) {
    231                 color = colorbar[colorbar.length-i-1];
    232                 color = [Math.round(color[0]*255),Math.round(color[1]*255),Math.round(color[2]*255)];   
    233                 tgradient.addColorStop(i/colorbar.length,'rgba('+color.toString()+',1.0)');
     228        if (options.getfieldvalue('log',10)!='off') {
     229                for (var i = 0; i < colorbar.length; i++) {
     230                        color = colorbar[colorbar.length-i-1];
     231                        color = [Math.round(color[0]*255),Math.round(color[1]*255),Math.round(color[2]*255)];   
     232                        tgradient.addColorStop(1.0-Math.exp(Math.log((colorbar.length-i)/colorbar.length)/Math.log(options.getfieldvalue('log',10))*i/colorbar.length),'rgba('+color.toString()+',1.0)');
     233                }
     234        } else {
     235                for (var i = 0; i < colorbar.length; i++) {
     236                        color = colorbar[colorbar.length-i-1];
     237                        color = [Math.round(color[0]*255),Math.round(color[1]*255),Math.round(color[2]*255)];   
     238                        tgradient.addColorStop(i/colorbar.length,'rgba('+color.toString()+',1.0)');
     239                }
    234240        }
    235241       

  • issm/trunk/src/m/plot/manualcb.m

    r20500 r20548  
    6262        ytick  = (ztick-zmin)/(zmax-zmin);
    6363else
    64         ztick = getfieldvalue(options,'tick',round(logspace(log(zmin)/log(10),log(zmax)/log(10),8)));
     64        %old method
     65        ztick = getfieldvalue(options,'tick',round( logspace(log10(zmin),log10(zmax),8) ));
    6566        ytick = linspace(0,1,numel(ztick));
     67
     68        %New method
     69        test=logspace(-10,10,21);
     70        pos=find(test>=zmin & test<=zmax);
     71        ztick= test(pos);
     72        ytick= (log(ztick) - log(zmin))/(log(zmax) - log(zmin));
    6673end
    6774
    … …  
    6976hold on
    7077numcolors=size(cmap,1);
    71 if strcmpi(getfieldvalue(options,'orientation','vertical'),'vertical'),
    72         image_rgb = ind2rgb(repmat((1:numcolors)',1,2),cmap);
    73 else
    74         image_rgb = ind2rgb(repmat((1:numcolors)',1,2)',cmap);
    75 end
    7678if 1,
    7779        %disappears somtimes
    78         imagesc([0 1],[0 1],repmat(image_rgb,[1 10 1]));
     80        if strcmpi(getfieldvalue(options,'orientation','vertical'),'vertical'),
     81                image_rgb = ind2rgb(repmat((1:numcolors)',1,10),cmap);
     82        else
     83                image_rgb = ind2rgb(repmat((1:numcolors),10,1),cmap);
     84        end
     85
     86        imagesc([0 1],[0 1],image_rgb);
    7987else
    8088        %Creates triangles when exported as pdf
    81         for i=1:numcolors,
    82                 patch([0,0,1,1],[(i-1)/numcolors,i/numcolors,i/numcolors,(i-1)/numcolors],'none','FaceColor',cmap(i,:),'Clipping','off','EdgeColor','none')
     89        if strcmpi(getfieldvalue(options,'orientation','vertical'),'vertical'),
     90                for i=1:numcolors,
     91                        patch([0,0,1,1],[(i-1)/numcolors,i/numcolors,i/numcolors,(i-1)/numcolors],0,'FaceColor',cmap(i,:),'Clipping','off','EdgeColor','none')
     92                end
     93        else
     94                for i=1:numcolors,
     95                        patch([(i-1)/numcolors,i/numcolors,i/numcolors,(i-1)/numcolors],[0,0,1,1],0,'FaceColor',cmap(i,:),'Clipping','off','EdgeColor','none')
     96                end
    8397        end
    8498end

  • issm/trunk/src/m/plot/plot_overlay.js

    r20500 r20548  
    127127                vertices[vertices.length] = vertex[2];
    128128               
    129                 texcoords[texcoords.length] = (x[i] - modelxlim[0]) / xrange;
    130                 texcoords[texcoords.length] = (y[i] - modelylim[0]) / yrange;
     129                if (md.mesh.classname() == 'mesh3dsurface') {
     130                        texcoords[texcoords.length] = degrees(Math.atan2(y[i], x[i])) / 360 + 0.5;
     131                        texcoords[texcoords.length] = degrees(Math.asin(z[i] / magnitude)) / 180 + 0.5;
     132                }
     133                else {
     134                        texcoords[texcoords.length] = (x[i] - modelxlim[0]) / xrange;
     135                        texcoords[texcoords.length] = (y[i] - modelylim[0]) / yrange;
     136                }
    131137        }
    132138

  • issm/trunk/src/m/plot/plot_unit.js

    r20500 r20548  
    9696                        else{ //triangular elements     
    9797                                caxis = options.getfieldvalue('caxis',[ArrayMin(data),ArrayMax(data)]);
     98                                if (options.getfieldvalue('log','off')!='off') caxis = [Math.log10(caxis[0])/Math.log10(options.getfieldvalue('log',10)),Math.log10(caxis[1])/Math.log10(options.getfieldvalue('log',10))];
    9899                                datamin = caxis[0];
    99100                                datamax = caxis[1];
    … …  
    196197                                        })
    197198                                });
    198                                
    199                                 for(var i = 0; i < data.length; i++){
    200                                         //Prevent evaluation of datasubarray min/max if caxis exists
    201                                         if (options.exist('caxis')) {
    202                                                 caxis = options.getfieldvalue('caxis');
    203                                         }
    204                                         else {
    205                                                 caxis = [ArrayMin(data[i]),ArrayMax(data[i].slice(0,-1))];
    206                                         }
     199                                //Prevent evaluation of datasubarray min/max if caxis exists
     200                                if (options.exist('caxis')) {
     201                                        caxis = options.getfieldvalue('caxis');
     202                                }
     203                                else {
     204                                        caxis = [ArrayMin(data[0]),ArrayMax(data[0].slice(0,-1))];
     205                                }
     206                                if (options.getfieldvalue('log','off')!='off') caxis = [Math.log10(caxis[0])/Math.log10(options.getfieldvalue('log',10)),Math.log10(caxis[1])/Math.log10(options.getfieldvalue('log',10))];
     207                               
     208                                for(var i = 0; i < data.length; i++){                                   
    207209                                        datamin = caxis[0];
    208210                                        datamax = caxis[1];

  • issm/trunk/src/m/plot/processdata.js

    r20500 r20548  
    143143
    144144                //log?
    145                 if (options.exist('log')){
     145                if (options.getfieldvalue('log','off')!='off'){
    146146                        var bounds=options.getfieldvalue('caxis',[ArrayMin(data),ArrayMax(data)]);
    147147                        for(var i=0;i<md.mesh.numberofelements;i++)if(data[i]<bounds[0])data[i]=bounds[0];
    … …  
    151151                        for(var i=0;i<md.mesh.numberofelements;i++){
    152152                                if(!IsNaN(data[i])){
    153                                         data[i]=Math.log10(data[i])/Math.log10(options.getfieldvalue('log'));
     153                                        data[i]=Math.log10(data[i])/Math.log10(options.getfieldvalue('log',10));
    154154                                }
    155155                        }
    … …  
    184184
    185185                //log?
    186                 if (options.exist('log')){
     186                if (options.getfieldvalue('log','off')!='off'){
    187187                        var bounds=options.getfieldvalue('caxis',[ArrayMin(data),ArrayMax(data)]);
    188188                        for(var i=0;i<md.mesh.numberofvertices;i++)if(data[i]<bounds[0])data[i]=bounds[0];
    … …  
    191191                                throw Error("Log option cannot be applied on negative values. Use caxis option (Rignot''s settings: [1.5 max(data)])");
    192192                        }
    193                    for(var i=0;i<md.mesh.numberofvertices;i++){
    194                            data[i]=Math.log10(data[i])/Math.log10(options.getfieldvalue('log'));
    195                    }
     193                        for(var i=0;i<md.mesh.numberofvertices;i++){
     194                           data[i]=Math.log10(data[i])/Math.log10(options.getfieldvalue('log',10));
     195                        }
    196196                }
    197197        }
    … …  
    200200    if (datasize==md.mesh.numberofvertices+1){
    201201        datatype=5;
     202               
     203                //log? 
     204                if (options.getfieldvalue('log','off')!='off'){
     205                        var bounds=options.getfieldvalue('caxis',[ArrayMin(data),ArrayMax(data)]);
     206                        for(var i=0;i<md.mesh.numberofvertices;i++) {
     207                                for(var j=0;j<data[i].length;j++) {
     208                                        if(data[i][j]<bounds[0])data[i][j]=bounds[0];
     209                                }
     210                        }
     211                        for(var i=0;i<md.mesh.numberofvertices;i++) {
     212                                for(var j=0;j<data[i].length;j++) {
     213                                        if(data[i][j]>bounds[1])data[i][j]=bounds[1];
     214                                }
     215                        }
     216                        for(var i=0;i<md.mesh.numberofvertices;i++) {
     217                                for(var j=0;j<data[i].length;j++) {
     218                                        if(data[i][j]<=0) {
     219                                                throw Error("Log option cannot be applied on negative values. Use caxis option (Rignot''s settings: [1.5 max(data)])");
     220                                        }
     221                                }
     222                        }
     223                        for(var i=0;i<md.mesh.numberofvertices;i++){
     224                                for(var j=0;j<data[i].length;j++) {
     225                                        data[i][j]=Math.log10(data[i][j])/Math.log10(options.getfieldvalue('log',10));
     226                                }
     227                        }
     228                }
    202229    }
    203230       

  • issm/trunk/src/wrappers/matlab/Makefile.am

    r20500 r20548  
    6161                                                 PropagateFlagsFromConnectivity.la\
    6262                                                 StringToEnum.la\
    63                                                  TriaSearch.la\
    6463                                                 TriMesh.la\
    6564                                                 TriMeshProcessRifts.la\
    … …  
    250249Shp2Kml_la_LIBADD = ${deps} $(SHAPELIBLIB) $(MPILIB) $(PETSCLIB) $(GSLLIB) $(PROJ4LIB)
    251250
    252 TriaSearch_la_SOURCES = ../TriaSearch/TriaSearch.cpp
    253 TriaSearch_la_LIBADD = ${deps} $(MPILIB) $(PETSCLIB) $(GSLLIB) $(PROJ4LIB)
    254 
    255251TriMesh_la_SOURCES = ../TriMesh/TriMesh.cpp
    256252TriMesh_la_LIBADD = ${deps} $(MPILIB) $(PETSCLIB) $(TRIANGLELIB) $(GSLLIB) $(PROJ4LIB)

  • issm/trunk/test

  • issm/trunk/test/NightlyRun/test124.m

    r20500 r20548  
    2121        'Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3'};
    2222field_tolerances={...
    23         2e-09,3.0e-9,3.0e-9,1e-10,1e-13,1e-12,1e-12,...
    24         1e-09,3e-10,8e-10,1e-10,1e-10,1e-10,1e-10,...
    25         1e-09,3e-10,8e-10,1e-10,1e-10,1e-10,1e-10};
     23        2e-09,3e-9,3.0e-9,1e-10,1e-13,1e-12,1e-12,...
     24        2e-09,3e-9,8e-10,1e-10,1e-10,1e-10,1e-10,...
     25        2e-09,3e-9,8e-10,1e-10,1e-10,1e-10,1e-10};
    2626field_values={...
    2727        (md.results.TransientSolution(1).Vx),...

  • issm/trunk/test/NightlyRun/test211.py

    r19105 r20548  
    2828        'Vx3','Vy3','Vz3','Vel3','Pressure3','Bed3','Surface3','Thickness3','Temperature3','BasalforcingsGroundediceMeltingRate3']
    2929field_tolerances=[\
    30         2e-08,2e-08,1e-06,2e-08,1e-08,1e-08,1e-08,1e-08,1e-08,1e-08,\
    31         5e-07,5e-07,5e-05,1e-07,1e-08,1e-08,1e-08,1e-08,5e-08,2e-06,\
    32         5e-06,5e-06,5e-05,5e-06,1e-07,5e-07,5e-07,5e-07,5e-06,5e-05]
     30        2e-08,2e-08,2e-06,2e-08,1e-08,1e-08,1e-08,1e-08,1e-08,1e-08,\
     31        5e-07,5e-07,5e-05,2e-07,1e-08,1e-08,1e-08,1e-08,6e-08,2e-06,\
     32        5e-06,5e-06,5e-05,5e-06,5e-07,5e-07,5e-07,5e-07,5e-06,5e-05]
    3333field_values=[\
    3434        md.results.TransientSolution[0].Vx,\

  • issm/trunk/test/NightlyRun/test3020.py

    r20500 r20548  
    6262
    6363md=solve(md,TransientSolutionEnum())
    64 V0=md.results.TransientSolution[2].IceVolume
    65 MaxV0=md.results.TransientSolution[2].MaxVel
     64V0=md.results.TransientSolution[-1].IceVolume
     65MaxV0=md.results.TransientSolution[-1].MaxVel
    6666
    6767#backward
    … …  
    7474
    7575md=solve(md,TransientSolutionEnum())
    76 V2=md.results.TransientSolution[2].IceVolume
    77 MaxV2=md.results.TransientSolution[2].MaxVel
     76V2=md.results.TransientSolution[-1].IceVolume
     77MaxV2=md.results.TransientSolution[-1].MaxVel
    7878
    7979#compute resulting derivative

  • issm/trunk/test/NightlyRun/test404.m

    r20500 r20548  
    1010%Fields and tolerances to track changes
    1111field_names     ={'Vx','Vy','Vz','Vel','Pressure'};
    12 field_tolerances={6e-07,6e-07,2e-06,1e-06,3e-07};
     12field_tolerances={6e-07,6e-07,2e-06,1e-06,5e-07};
    1313field_values={...
    1414        (md.results.StressbalanceSolution.Vx),...

  • issm/trunk/test/NightlyRun/test404.py

    r20500 r20548  
    2020#Fields and tolerances to track changes
    2121field_names     =['Vx','Vy','Vz','Vel','Pressure']
    22 field_tolerances=[6e-07,6e-07,2e-06,1e-06,3e-07]
     22field_tolerances=[6e-07,6e-07,2e-06,1e-06,5e-07]
    2323field_values=[\
    2424        md.results.StressbalanceSolution.Vx,\

  • issm/trunk/test/NightlyRun/test422.m

    r20500 r20548  
    1111%Fields and tolerances to track changes
    1212field_names     ={'Vx','Vy','Vz','Vel','Pressure'};
    13 field_tolerances={4e-07,4e-07,2e-06,2e-07,5e-07};
     13field_tolerances={4e-07,4e-07,2e-06,4e-07,5e-07};
    1414field_values={...
    1515        (md.results.StressbalanceSolution.Vx),...

  • issm/trunk/test/NightlyRun/test422.py

    r20500 r20548  
    2121#Fields and tolerances to track changes
    2222field_names     =['Vx','Vy','Vz','Vel','Pressure']
    23 field_tolerances=[4e-07,4e-07,2e-06,2e-07,5e-07]
     23field_tolerances=[4e-07,4e-07,2e-06,4e-07,5e-07]
    2424field_values=[\
    2525        md.results.StressbalanceSolution.Vx,\
Note: See TracChangeset for help on using the changeset viewer.