Changeset 25836 for issm/trunk/src/m/classes/sealevelmodel.m

Timestamp:

12/08/20 08:45:53 (4 years ago)

Author:

Mathieu Morlighem

Message:

merged trunk-jpl and trunk for revision 25834

Location:

issm/trunk

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/m/classes/sealevelmodel.m (modified) (15 diffs)

issm/trunk/src/m/classes/sealevelmodel.m

@@ -12 +12 @@
         properties (SetAccess=public) %Model fields
                 % {{{
-                icecaps          = {}; % list of ice cap models
+                icecaps          = {}; % list of land/ice models, name should  be change longer term.
                 earth            = 0;  % model for the whole earth
                 basins           = {}; % list  of basins, matching icecaps, where shapefile info is held.
@@ -19 +19 @@
                 settings         = 0;
                 private          = 0;
-                range            = 0;
                 mergedcaps       = 0;
                 transitions      = {};
-                eltransitions      = {};
+                eltransitions    = {};
+                planet           = '';
                 %}}}
         end
         methods
                 function slm = sealevelmodel(varargin) % {{{
-
-                        if nargin==0, 
-                                slm=setdefaultparameters(slm);
-                        else 
-                                slm=setdefaultparameters(slm);
-
-                                options=pairoptions(varargin{:}); 
-                        
+                        slm=setdefaultparameters(slm);
+
+                        if nargin==1,
+
+                                options=pairoptions(varargin{:});
+
                                 %recover all the icecap models: 
                                 slm.icecaps=getfieldvalues(options,'ice_cap',{}); 
                                 
                                 %recover the earth model:
-                                slm.earth = getfieldvalue(options,'earth');
+                                slm.earth = getfieldvalue(options,'earth',0);
+
+                                %set planet type: 
+                                slm.planet=getfieldvalue(options,'planet','earth');
+
                         end
                 end
@@ -89 +91 @@
                         slm.private           = private();
                         slm.cluster           = generic();
-                        slm.range             = {};
                         slm.transitions       = {};
-                        slm.eltransitions       = {};
+                        slm.eltransitions     = {};
+                        slm.planet            = 'earth';
                 end
                 %}}}
@@ -100 +102 @@
                         disp(sprintf('%19s: %-22s -- %s','cluster'         ,['[1x1 ' class(self.cluster) ']'],'cluster parameters (number of cpus...)'));
                         disp(sprintf('%19s: %-22s -- %s','miscellaneous'   ,['[1x1 ' class(self.miscellaneous) ']'],'miscellaneous fields'));
-                        disp(sprintf('%19s: %-22s -- %s','range'   ,['[1x1 ' class(self.range) ']'],'ranges'));
                 end % }}}
                 function self=mergeresults(self) % {{{
-                        champs=fields(self.icecaps{1}.results.TransientSolution);
+                        champs=fieldnames(self.icecaps{1}.results.TransientSolution);
                         for i=1:length(self.mergedcaps)/2,
                                 md=self.mergedcaps{2*(i-1)+1}; trans=self.mergedcaps{2*(i-1)+2};
-                                icecaps=self.icecaps(self.range{2*(i-1)+2});
+                                %icecaps=self.icecaps(self.range{2*(i-1)+2});
                                 for j=1:length(self.icecaps{1}.results.TransientSolution),
                                         for k=1:length(champs),
@@ -143 +144 @@
                         end
                 end % }}}
+                function n=ncaps(self) % {{{
+                        n=length(self.icecaps);
+                end % }}}
+                function list=continents(self) % {{{
+                        list={};
+                        for  i=1:length(self.basins),
+                                list{end+1}=self.basins{i}.continent;
+                        end
+                        list=unique(list);
+                end % }}}
+                function list=basinsfromcontinent(self,continent) % {{{
+                        list={};
+                        for  i=1:length(self.icecaps),
+                                if strcmpi(self.basins{i}.continent,continent),
+                                        list{end+1}=self.basins{i}.name;
+                                end
+                        end
+                        list=unique(list);
+                end % }}}
                 function addbasin(self,bas) % {{{
-                if ~strcmpi(class(bas),'basin')
-                        error('addbasin method only takes a ''basin'' class object as input');
-                end;
-                self.basins{end+1}=bas;
+                        if ~strcmpi(class(bas),'basin')
+                                error('addbasin method only takes a ''basin'' class object as input');
+                        end;
+                        self.basins{end+1}=bas;
+                end % }}}
+                function intersections2d(self,varargin) % {{{
+
+                        options=pairoptions(varargin{:}); 
+                        force=getfieldvalue(options,'force',0);
+                        
+                        %initialize, to avoid issues of having more transitions than meshes. 
+                        self.transitions={};
+                        self.eltransitions={};
+
+                        %for elements:
+                        xe=self.earth.mesh.x(self.earth.mesh.elements)*[1;1;1]/3;
+                        ye=self.earth.mesh.y(self.earth.mesh.elements)*[1;1;1]/3;
+                        
+                        for i=1:length(self.icecaps),
+                                mdi=self.icecaps{i};
+                
+                                %for elements: 
+                                xei=mdi.mesh.x(mdi.mesh.elements)*[1;1;1]/3;
+                                yei=mdi.mesh.y(mdi.mesh.elements)*[1;1;1]/3;
+                
+                                disp(sprintf('Computing vertex intersections for basin %s',self.basins{i}.name));
+                        
+                                self.transitions{end+1}=meshintersect2d(self.earth.mesh.x,self.earth.mesh.y,mdi.mesh.x,mdi.mesh.y,'force',force);
+
+                                self.eltransitions{end+1}=meshintersect2d(xe,ye,xei,yei,'force',force);
+                        end
                 end % }}}
                 function intersections(self,varargin) % {{{
 
-                options=pairoptions(varargin{:}); 
-                force=getfieldvalue(options,'force',0);
-
-                %for elements:
-                xe=self.earth.mesh.x(self.earth.mesh.elements)*[1;1;1]/3;
-                ye=self.earth.mesh.y(self.earth.mesh.elements)*[1;1;1]/3;
-                ze=self.earth.mesh.z(self.earth.mesh.elements)*[1;1;1]/3;
-
-                for i=1:length(self.icecaps),
-                        mdi=self.icecaps{i};
-                        mdi=TwoDToThreeD(mdi);
-        
-                        %for elements: 
-                        xei=mdi.mesh.x(mdi.mesh.elements)*[1;1;1]/3;
-                        yei=mdi.mesh.y(mdi.mesh.elements)*[1;1;1]/3;
-                        zei=mdi.mesh.z(mdi.mesh.elements)*[1;1;1]/3;
-        
-                        disp(sprintf('Computing vertex intersections for basin %s',self.basins{i}.name));
+                        options=pairoptions(varargin{:}); 
+                        force=getfieldvalue(options,'force',0);
+                        
+                        %initialize, to avoid issues of having more transitions than meshes. 
+                        self.transitions={};
+                        self.eltransitions={};
+
+                        %for elements:
+                        xe=self.earth.mesh.x(self.earth.mesh.elements)*[1;1;1]/3;
+                        ye=self.earth.mesh.y(self.earth.mesh.elements)*[1;1;1]/3;
+                        ze=self.earth.mesh.z(self.earth.mesh.elements)*[1;1;1]/3;
+                        
+                        for i=1:length(self.icecaps),
+                                mdi=self.icecaps{i};
+                                mdi=TwoDToThreeD(mdi,self.planet);
                 
-                        self.transitions{end+1}=meshintersect3d(self.earth.mesh.x,self.earth.mesh.y,self.earth.mesh.z,mdi.mesh.x,mdi.mesh.y,mdi.mesh.z,'force',force);
-
-                        self.eltransitions{end+1}=meshintersect3d(xe,ye,ze,xei,yei,zei,'force',force);
-                end
-
+                                %for elements: 
+                                xei=mdi.mesh.x(mdi.mesh.elements)*[1;1;1]/3;
+                                yei=mdi.mesh.y(mdi.mesh.elements)*[1;1;1]/3;
+                                zei=mdi.mesh.z(mdi.mesh.elements)*[1;1;1]/3;
+                
+                                disp(sprintf('Computing vertex intersections for basin %s',self.basins{i}.name));
+                        
+                                self.transitions{end+1}=meshintersect3d(self.earth.mesh.x,self.earth.mesh.y,self.earth.mesh.z,mdi.mesh.x,mdi.mesh.y,mdi.mesh.z,'force',force);
+
+                                self.eltransitions{end+1}=meshintersect3d(xe,ye,ze,xei,yei,zei,'force',force);
+                        end
                 end % }}}
                 function checkintersections(self) % {{{
-                flags=zeros(self.earth.mesh.numberofvertices,1);
-                for i=1:length(self.basins),
-                        flags(self.transitions{i})=i;
-                end
-                plotmodel(self.earth,'data',flags,'coastline','on');
+                        flags=zeros(self.earth.mesh.numberofvertices,1);
+                        for i=1:length(self.basins),
+                                flags(self.transitions{i})=i;
+                        end
+                        plotmodel(self.earth,'data',flags,'coastline','on');
+
+                end % }}}
+                function checkbasinconsistency(self) % {{{
+                        for i=1:length(self.basins),
+                                self.basins{i}.checkconsistency();
+                        end
 
                 end % }}}
@@ -192 +248 @@
                         if iscell(continent),
                                 if length(continent)==1,
-                                         if strcmpi(continent{1},'all'),
-                                                 %need to transform this into a list of continents: 
-                                                 continent={}; 
-                                                 for i=1:length(self.basins),
-                                                         continent{end+1}=self.basins{i}.continent;
-                                                 end
-                                                 continent=unique(continent);
-                                         end
+                                        if strcmpi(continent{1},'all'),
+                                                %need to transform this into a list of continents: 
+                                                continent={}; 
+                                                for i=1:length(self.basins),
+                                                        continent{end+1}=self.basins{i}.continent;
+                                                end
+                                                continent=unique(continent);
+                                        else
+                                                %nothing to do: assume we have a list of continents
+                                        end
                                 else
-                                        %nothing to do, we have a list of continents
+                                        %nothing to do: assume we have a list of continents
                                 end
                         else
@@ -208 +266 @@
                                         continent={}; 
                                         for i=1:length(self.basins),
-                                                 continent{end+1}=self.basins{i}.continent;
+                                                continent{end+1}=self.basins{i}.continent;
                                         end
                                         continent=unique(continent);
@@ -219 +277 @@
                         if iscell(bas),
                                 if length(bas)==1,
-                                         if strcmpi(bas{1},'all'),
-                                                 %need to transform this into a list of basins: 
-                                                 baslist=[];
-                                                 for i=1:length(self.basins),
-                                                         if self.basins{i}.iscontinentany(continent{:}),
-                                                                 baslist(end+1)=i;
-                                                         end
-                                                 end
-                                                 baslist=unique(baslist);
-                                         else
-                                                 bas=bas{1};
-                                                 baslist=[];
-                                                 for i=1:length(self.basins),
-                                                         if self.basins{i}.iscontinentany(continent{:}),
-                                                                 if self.basins{i}.isnameany(bas),
-                                                                         baslist(end+1)=i;
-                                                                 end
-                                                         end
-                                                 end
-
-                                         end
+                                        if strcmpi(bas{1},'all'),
+                                                %need to transform this into a list of basins: 
+                                                baslist=[];
+                                                for i=1:length(self.basins),
+                                                        if self.basins{i}.iscontinentany(continent{:}),
+                                                                baslist(end+1)=i;
+                                                        end
+                                                end
+                                                baslist=unique(baslist);
+                                        else
+                                        bas=bas{1};
+                                        baslist=[];
+                                        for i=1:length(self.basins),
+                                                if self.basins{i}.iscontinentany(continent{:}),
+                                                        if self.basins{i}.isnameany(bas),
+                                                                baslist(end+1)=i;
+                                                        end
+                                                end
+                                        end
+
+                                        end
                                 else
                                         %we have a list of basin names: 
@@ -280 +338 @@
                 end % }}}
                 function addicecap(self,md) % {{{
-                if ~strcmpi(class(md),'model')
-                        error('addicecap method only takes a ''model'' class object as input');
-                end
-                self.icecaps{end+1}=md;
+                        if ~strcmpi(class(md),'model')
+                                error('addicecap method only takes a ''model'' class object as input');
+                        end
+                        self.icecaps{end+1}=md;
                 end % }}}
                 function basinsplot3d(self,varargin) % {{{
-                for i=1:length(self.basins),
-                        self.basins{i}.plot3d(varargin{:});
-                end
+                        for i=1:length(self.basins),
+                                self.basins{i}.plot3d(varargin{:});
+                        end
                 end % }}}
                 function caticecaps(self,varargin) % {{{
@@ -300 +358 @@
                         models=self.icecaps;
                         for i=1:length(models),
-                                models{i}=TwoDToThreeD(models{i});
+                                models{i}=TwoDToThreeD(models{i},self.planet);
                         end
                         
@@ -326 +384 @@
                         self.earth=md;
 
+                        %Create mesh radius:
+                        self.earth.mesh.r=planetradius('earth')*ones(md.mesh.numberofvertices,1);
+
+                end % }}}
+                function caticecaps2d(self,varargin) % {{{
+
+                        %recover options: 
+                        options=pairoptions(varargin{:}); 
+                        tolerance=getfieldvalue(options,'tolerance',1e-5);
+                        loneedgesdetect=getfieldvalue(options,'loneedgesdetect',0);
+                        models=self.icecaps;
+
+                        %Plug all models together:
+                        md=models{1}; 
+                        for i=2:length(models),
+                                md=modelmerge2d(md,models{i},'tolerance',tolerance);
+                        end
+
+                        %Look for lone edges if asked for it: {{{
+                        if loneedgesdetect,
+                                edges=loneedges(md);
+                                hold on;
+                                for i=1:length(edges),
+                                        ind1=edges(i,1);
+                                        ind2=edges(i,2);
+                                        plot([md.mesh.x(ind1),md.mesh.x(ind2)],[md.mesh.y(ind1),md.mesh.y(ind2)],'g*-');
+                                end
+                        end %}}}
+
+                        %Plug into earth: 
+                        self.earth=md;
+
                 end % }}}
                 function viscousiterations(self) % {{{
                         for  i=1:length(self.icecaps),
-                                ic=self.icecaps{i}; 
+                                ic=self.icecaps{i};
                                 mvi=ic.results.TransientSolution(1).StressbalanceConvergenceNumSteps;
                                 for j=2:length(ic.results.TransientSolution)-1,
@@ -339 +429 @@
                 function maxtimestep(self) % {{{
                         for  i=1:length(self.icecaps),
-                                ic=self.icecaps{i}; 
+                                ic=self.icecaps{i};
                                 mvi=length(ic.results.TransientSolution);
                                 timei=ic.results.TransientSolution(end).time;
@@ -350 +440 @@
                 function transfer(self,string) % {{{
                         %Recover field size in one icecap: 
-                        eval(['n=length(self.icecaps{1}.' string ');']);
+                        eval(['n=size(self.icecaps{1}.' string ',1);']);
                         if n==self.icecaps{1}.mesh.numberofvertices,
                                 eval(['self.earth.' string '=zeros(self.earth.mesh.numberofvertices,1);']);
@@ -356 +446 @@
                                         eval(['self.earth.' string '(self.transitions{' num2str(i) '})=self.icecaps{' num2str(i) '}.' string ';']);
                                 end
-                        elseif n==self.icecaps{1} .mesh.numberofelements,
+                        elseif n==(self.icecaps{1}.mesh.numberofvertices+1),
+                                %dealing with a transient dataset. 
+                                %check that all timetags are similar between all icecaps:  %{{{
+                                for i=1:length(self.icecaps),
+                                        eval(['capfieldi= self.icecaps{' num2str(i) '}.' string ';']);
+                                        for j=(i+1):length(self.icecaps),
+                                                eval(['capfieldj= self.icecaps{' num2str(j) '}.' string ';']);
+                                                if ~isequal(capfieldi(end,:),capfieldj(end,:)),
+                                                        error(['Time stamps for ' string ' field are different between icecaps ' num2str(i) ' and ' num2str(j)]);
+                                                end
+                                        end
+                                end
+                                eval(['capfield1= self.icecaps{1}.' string ';']); 
+                                times=capfield1(end,:);
+                                nsteps=length(times);
+                                %}}}
+                                %initialize:  %{{{
+                                eval(['field=zeros(self.earth.mesh.numberofvertices+1,' num2str(nsteps) ');']);
+                                eval(['field(end,:)=times;']); %transfer the times only, not the values
+                                %}}}
+                                %transfer all time fields: {{{
+                                for i=1:length(self.icecaps),
+                                        eval(['capfieldi= self.icecaps{' num2str(i) '}.' string ';']);
+                                        for j=1:nsteps,
+                                                eval(['field(self.transitions{' num2str(i) '},' num2str(j) ')=capfieldi(1:end-1,' num2str(j) ');']); %transfer only the values, not the time.
+                                        end
+                                end
+                                eval(['self.earth.' string '=field;']); %do not forget to plug the field variable into its final location
+                                %}}}
+                        elseif n==self.icecaps{1}.mesh.numberofelements,
                                 eval(['self.earth.' string '=zeros(self.earth.mesh.numberofelements,1);']);
                                 for i=1:length(self.icecaps),
@@ -389 +508 @@
                         self.earth=ic;
                 end % }}}
+                function self=initializemodels(self) % {{{
+
+                        for i=1:length(self.basins),
+                                md=model();
+                                md.miscellaneous.name=self.basins{i}.name;
+                                self.addicecap(md);
+                        end
+                end % }}}
         end
 end