Index: /issm/trunk-jpl/src/m/classes/modellist.m =================================================================== --- /issm/trunk-jpl/src/m/classes/modellist.m (revision 13007) +++ /issm/trunk-jpl/src/m/classes/modellist.m (revision 13008) @@ -10,4 +10,116 @@ end methods + function md_list=modelsextract(md,flags,minel,varargin) % {{{ + %modelsextract - extract several self contained models according to a list of element flags. + % + % The difference between this routine and the modelextract.m routine (without an 's') is that + % as many models are extracted as there are closed contours defined in area. + % This routine is needed for example when doing data assimilation of ice shelves in Antarctica. + % Many independent ice shelves are present, and we don't want data assimilation on one ice shelf + % to be hindered by another totally independent ice shelf. + % + % Usage: + % md_list=modelsextract(md,elementfalgs,minel); + % + % Examples: + % md_list=modelsextract(md,md.mask.elementonfloatingice,1000); + % + % See also: EXTRUDE, COLLAPSE, MODELEXTRACT + + disp('selecting pools of elements'); + %go through flags and build as many independent element flags as there are groups of connected 1s + %in flags. + + %2D or 3D? + if md.mesh.dimension==3, + numberofelements=md.mesh.numberofelements2d; %this will be forgotten when we get out. + flags=project2d(md,flags,1); + else + numberofelements=md.mesh.numberofelements; + end + + %recover extra arguments: + distance=0; + if nargin==4, + distance=varargin{1}; + end + + flag_list=cell(0,1); + + for i=1:size(flags,1), + + if (flags(i)), + + %ok, we are sure element i is part of a new pool. + pool=zeros(numberofelements,1); + pool=PropagateFlagsFromConnectivity(md.mesh.elementconnectivity,pool,i,flags); + flag_list{end+1,1}=pool; + + %speed up rest of computation by taking pool out of flags: + pos=find(pool);flags(pos)=0; + + end + end + + %go through flag_list and discard any pool of less than minel elements: + ex_pos=[]; + for i=1:length(flag_list), + if length(find(flag_list{i}))=1)); + +%average x and y over elements: +x_elem=md.mesh.x(md.mesh.elements)*[1;1;1]/3; +y_elem=md.mesh.y(md.mesh.elements)*[1;1;1]/3; + +while 1, + + %keep copy of new_flags for this loop: + new_flags_bak=new_flags; + + %extend new flags by connectivity + pos=find(new_flags); + + connected_elements=md.mesh.elementconnectivity(pos,:); + connected_elements=connected_elements(find(connected_elements)); + new_flags(connected_elements)=1; + + %get new elements: + new_elements=find(new_flags & ~new_flags_bak); + if ~length(new_elements), + %we are done! + break; + end + + %check which of these new elements are more than distance away from the border elements + for i=1:length(new_elements), + dist=sqrt( (x_elem(border_elements)-x_elem(new_elements(i))).^2 + (y_elem(border_elements)-y_elem(new_elements(i))).^2)-distance; + if ~any(dist<0) + %none of the border elements are within distance, this element is outside out area of interest. + %ensure this element never gets found again in the connectivity. + pos=find(md.mesh.elementconnectivity==new_elements(i)); + md.mesh.elementconnectivity(pos)=0; + %exclude this new element from the new_flags! + new_flags(new_elements(i))=0; + end + end +end Index: /issm/trunk-jpl/src/m/contrib/hack/thicknessevolution.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/hack/thicknessevolution.m (revision 13008) +++ /issm/trunk-jpl/src/m/contrib/hack/thicknessevolution.m (revision 13008) @@ -0,0 +1,28 @@ +function dhdt=thicknessevolution(md) +%THICKNESSEVOLUTION - compute the new thickness of a model after ∆t +% +% This routine compute the new thickness of a model after a time step +% according to the following formula: +% dh/dt=-div(Hu) +% +% Usage: +% dhdt=thicknessevolution(md) + +if (length(md.initialization.vx)~=md.mesh.numberofvertices)|(length(md.initialization.vy)~=md.mesh.numberofvertices) + error(['thicknessevolution error message: vx and vy should have a length of ' num2str(md.mesh.numberofvertices)]) +end + +%load some variables +H=md.geometry.thickness; +vx=md.initialization.vx; +vy=md.initialization.vy; +index=md.mesh.elements; + +%compute nodal functions coefficients N(x,y)=alpha x + beta y + gamma +[alpha beta]=GetNodalFunctionsCoeff(md.mesh.elements,md.mesh.x,md.mesh.y); + +%compute dhdt=div(Hu) +summation=1/3*ones(3,1); +dhdt=(vx(index)*summation).*sum( H(index).*alpha,2) + (vy(index)*summation).*sum(H(index).*beta,2) ... + + ( H(index)*summation).*sum(vx(index).*alpha,2) + ( H(index)*summation).*sum(vy(index).*beta,2); +dhdt=-dhdt; Index: /issm/trunk-jpl/src/m/contrib/hack/tres.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/hack/tres.m (revision 13008) +++ /issm/trunk-jpl/src/m/contrib/hack/tres.m (revision 13008) @@ -0,0 +1,119 @@ +function md=tres(md,string) +%TRES - transfer results results to corresponding model fields. +% +% Usage: md=tres(md,string) +% +% Example: md=tres(md,'diagnostic'); + +%check number of arguments + +if strcmpi(string,'diagnostic'), + if md.mesh.dimension==2, + md.initialization.vx=md.results.DiagnosticSolution.Vx; + md.initialization.vy=md.results.DiagnosticSolution.Vy; + else + md.initialization.vx=md.results.DiagnosticSolution.Vx; + md.initialization.vy=md.results.DiagnosticSolution.Vy; + md.initialization.vz=md.results.DiagnosticSolution.Vz; + end + md.initialization.vel=md.results.DiagnosticSolution.Vel; + + if isfield(md.results.DiagnosticSolution,'Pressure'), + md.initialization.pressure=md.results.DiagnosticSolution.Pressure; + end + if md.rifts.numrifts, + if isfield(md.results.DiagnosticSolution,'riftproperties'), + md.rifts.riftproperties=md.results.DiagnosticSolution.riftproperties; + end + end + if md.inversion.iscontrol==1, + for control_parameters=md.inversion.control_parameters + %Will need to be updated... good luck ;) + md.(EnumToModelField(control_parameters))=md.results.DiagnosticSolution.(EnumToString(control_parameters)); + end + end + +elseif strcmpi(string,'dakota'), + md.qmu.results=md.results.dakota; + +elseif strcmpi(string,'flaim'), + md.flaim.solution=md.results.FlaimSolution.solution; + md.flaim.quality =md.results.FlaimSolution.quality; + +elseif strcmpi(string,'transient'), + results=md.results.TransientSolution; + results2.Vel=NaN; + count=1; + for i=1:length(results), + if ~isempty(md.results.TransientSolution(i).Vel), + results2(count).Vel=md.results.TransientSolution(i).Vel; + results2(count).Surface=md.results.TransientSolution(i).Surface; + results2(count).Thickness=md.results.TransientSolution(i).Thickness; + results2(count).Bed=md.results.TransientSolution(i).Bed; + results2(count).Vx=md.results.TransientSolution(i).Vx; + results2(count).Vy=md.results.TransientSolution(i).Vy; + results2(count).time=md.results.TransientSolution(i).time; + results2(count).step=md.results.TransientSolution(i).step; + if ~strcmpi(md.groundingline.migration,'None'), + results2(count).ElementOnIceShelf=md.results.TransientSolution(i).ElementOnIceShelf; + end + count=count+1; + end + end + md.results.TransientSolution=results2; + clear results,results2; +elseif strcmpi(string,'steadystate'), + md.initialization.vx=md.results.SteadystateSolution.Vx; + md.initialization.vy=md.results.SteadystateSolution.Vy; + if isfield(md.results.SteadystateSolution,'Vz'), + md.initialization.vz=md.results.SteadystateSolution.Vz; + end + + md.initialization.vel=md.results.SteadystateSolution.Vel; + md.initialization.pressure=md.results.SteadystateSolution.Pressure; + md.initialization.temperature=md.results.SteadystateSolution.Temperature; + md.basalforcings.melting_rate=md.results.SteadystateSolution.BasalforcingsMeltingRate; + + if md.inversion.iscontrol==1, + for control_parameters=md.inversion.control_parameters + md.(EnumToModelField(control_parameters))=md.results.SteadystateSolution.(EnumToString(control_parameters)); + end + end + +elseif strcmpi(string,'thermal'), + md.initialization.temperature=md.results.ThermalSolution.Temperature; + md.basalforcings.melting_rate=md.results.ThermalSolution.BasalMeltingRate; +elseif strcmpi(string,'hydrology'), + md.initialization.watercolumn=md.results.HydrologySolution.Watercolumn; + +else + error(['tres error message: analysis ' string ' not supported yet!']); +end +end +function string=EnumToModelField(enum) % {{{ + %ENUMTOMODELFIELD - output string of model field associated to enum + % + % Usage: + % string=EnumToModelField(enum) + + disp('Warning: EnumToModelField is deprecated, it cannot work with new model definition. This function will be removed in the future'); + + switch enum, + + case ThicknessEnum(), string='thickness'; return + case FrictionCoefficientEnum(), string='drag_coefficient'; return + case MaterialsRheologyBEnum(), string='rheology_B'; return + case MaterialsRheologyBbarEnum(), string='rheology_B'; return + case BalancethicknessThickeningRateEnum: string='dhdt'; return + case VxEnum(), string='vx'; return + case InversionVxObsEnum(), string='vx_obs'; return + case VyEnum(), string='vy'; return + case InversionVyObsEnum(), string='vy_obs'; return + case BasalforcingsMeltingRateEnum(), string='basal_melting_rate'; return + case SurfaceforcingsAccumulationRateEnum(), string='surface_accumulation_rate'; return + case SurfaceforcingsAblationRateEnum(), string='surface_ablation_rate'; return + case SurfaceforcingsMassBalanceEnum(), string='surface_mass_balance'; return + otherwise, error(['Enum ' num2str(enum) ' not found associated to any model field']); + + end + end % }}} Index: /issm/trunk-jpl/src/m/contrib/massbalance/outflow.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/massbalance/outflow.m (revision 13008) +++ /issm/trunk-jpl/src/m/contrib/massbalance/outflow.m (revision 13008) @@ -0,0 +1,18 @@ +function flag=outflow(md); +%OUTFLOW - flag nodes on outflux boundary +% +% Usage: +% flag=outflow(md); + +A=md.mesh.segments(:,1); +B=md.mesh.segments(:,2); +Nx=-(md.mesh.y(A)-md.mesh.y(B)); +Ny= md.mesh.x(A)-md.mesh.x(B); +Vx=(md.initialization.vx(A)+md.initialization.vx(B))/2; +Vy=(md.initialization.vy(A)+md.initialization.vy(B))/2; + +%dot product +VdotN=Vx.*Nx+Vy.*Ny; + +flag=zeros(md.mesh.numberofvertices,1); +flag(A(find(VdotN>0)))=1; Index: /issm/trunk-jpl/src/m/geometry/ThicknessCorrection.m =================================================================== --- /issm/trunk-jpl/src/m/geometry/ThicknessCorrection.m (revision 13008) +++ /issm/trunk-jpl/src/m/geometry/ThicknessCorrection.m (revision 13008) @@ -0,0 +1,77 @@ +function md=ThicknessCorrection(md,varargin) +%THICKNESSCORRECTION - correct the thickness of the ice shelf near the grounding line +% +% This routine corrects the thickness and the bed on the transition zone +% by forcing the hydrostatic equilibrium. +% the thickness is modified as follows: +% thickness = (1-coeff) * thickness_observation + coeff * thickness_hydrostatic +% where: +% coeff=(d/distance)^2; +% distance=10km by default but can be specified +% +% Usage: +% md=ThicknessCorrection(md,varargin); +% +% Example: +% md=ThicknessCorrection(md); +% md=ThicknessCorrection(md,15000); + +%initialize thickness with the observations, and get hydrostatic thickness from the dem +thickness=md.geometry.thickness; +thickness_hydro=md.geometry.surface/(1-md.materials.rho_ice/md.materials.rho_water); +hydrostatic_ratio=zeros(size(md.geometry.thickness)); + +%get nodes on ice sheet and on ice shelf +pos_shelf=find(~md.mask.vertexongroundedice); +pos_GL=intersect(unique(md.mesh.elements(find(md.mask.elementongroundedice),:)),unique(md.mesh.elements(find(md.mask.elementonfloatingice),:))); +debug=(length(pos_shelf)>50000); + +%check that there is a GL +if isempty(pos_GL) + error('ThicknessCorrection error message: no grounding line has been detected. Check the model mask'); +end + +%get distance +if nargin==2, + distance=varargin{1}; +else + distance=10000; +end + +%modify thickness +if (debug), fprintf('%s',' correction progress: 0.00 %'); end +for i=1:length(pos_shelf) + + if (debug & mod(i,100)==0), + fprintf('\b\b\b\b\b\b\b%5.2f%s',i/length(pos_shelf)*100,' %'); + end + + %search the node on ice sheet the closest to i + [d posd]=min(sqrt((md.mesh.x(pos_shelf(i))-md.mesh.x(pos_GL)).^2+(md.mesh.y(pos_shelf(i))-md.mesh.y(pos_GL)).^2)); + + if d>distance, + + %if d > 15km, hydrostatic equilibrium + hydrostatic_ratio(pos_shelf(i))=1; + thickness(pos_shelf(i))=thickness_hydro(pos_shelf(i)); + + else + + %else: quadratic combination of hydrostatic equilibrium and observations + hydrostatic_ratio(pos_shelf(i))=(d/distance)^2; + thickness(pos_shelf(i))=(1-hydrostatic_ratio(pos_shelf(i)))*thickness(pos_shelf(i))+hydrostatic_ratio(pos_shelf(i))*thickness_hydro(pos_shelf(i)); + + end +end +if (debug), fprintf('\b\b\b\b\b\b\b%5.2f%s\n',100,' %'); end + +%check the computed thickness +minth=1/(1-md.materials.rho_ice/md.materials.rho_water); +pos=find(isnan(thickness) | (thickness<=0)); +thickness(pos)=minth; +hydrostatic_ratio(pos)=-1; + +%change bed to take into account the changes in thickness +md.geometry.thickness=thickness; +md.geometry.hydrostatic_ratio=hydrostatic_ratio; +md.geometry.bed=md.geometry.surface-md.geometry.thickness; Index: /issm/trunk-jpl/src/m/geometry/recover_areas.m =================================================================== --- /issm/trunk-jpl/src/m/geometry/recover_areas.m (revision 13008) +++ /issm/trunk-jpl/src/m/geometry/recover_areas.m (revision 13008) @@ -0,0 +1,22 @@ +function [hutterflag macayealflag pattynflag stokesflag filltype]=recover_areas(md,varargin); +%RECOVER_AREAS - flag the element depending on the physical model that is assigned to them +% +% This routine is called by setelementstype, do not use +% +% Usage: +% [hutterflag macayealflag pattynflag stokesflag filltype]=recover_areas(md,varargin); + + %go through varargin, extract options and plug them into subtype options, by order of appearance + options=pairoptions(varargin{:}); + options=deleteduplicates(options,1); + + %recover elements distribution + hutterflag =FlagElements(md,getfieldvalue(options,'hutter','')); + macayealflag=FlagElements(md,getfieldvalue(options,'macayeal','')); + pattynflag =FlagElements(md,getfieldvalue(options,'pattyn','')); + stokesflag =FlagElements(md,getfieldvalue(options,'stokes','')); + filltype =getfieldvalue(options,'fill','none'); + +end %end function + + Index: /issm/trunk-jpl/src/m/geometry/recover_areas.py =================================================================== --- /issm/trunk-jpl/src/m/geometry/recover_areas.py (revision 13008) +++ /issm/trunk-jpl/src/m/geometry/recover_areas.py (revision 13008) @@ -0,0 +1,26 @@ +from pairoptions import * +from FlagElements import * + +def recover_areas(md,*args): + """ + RECOVER_AREAS - flag the element depending on the physical model that is assigned to them + + This routine is called by setelementstype, do not use + + Usage: + [hutterflag macayealflag pattynflag stokesflag filltype]=recover_areas(md,varargin); + """ + + #go through varargin, extract options and plug them into subtype options, by order of appearance + options=pairoptions(*args) +# options=deleteduplicates(options,1); + + #recover elements distribution + hutterflag =FlagElements(md,options.getfieldvalue('hutter','')) + macayealflag=FlagElements(md,options.getfieldvalue('macayeal','')) + pattynflag =FlagElements(md,options.getfieldvalue('pattyn','')) + stokesflag =FlagElements(md,options.getfieldvalue('stokes','')) + filltype =options.getfieldvalue('fill','none') + + return hutterflag,macayealflag,pattynflag,stokesflag,filltype + Index: /issm/trunk-jpl/src/m/geometry/slope.m =================================================================== --- /issm/trunk-jpl/src/m/geometry/slope.m (revision 13008) +++ /issm/trunk-jpl/src/m/geometry/slope.m (revision 13008) @@ -0,0 +1,32 @@ +function [sx,sy,s]=slope(md) +%SLOPE - compute the surface slope +% +% Usage: +% [sx,sy,s]=slope(md) + +%load some variables (it is much faster if the variab;es are loaded from md once for all) +if (md.mesh.dimension==2), + numberofelements=md.mesh.numberofelements; + numberofnodes=md.mesh.numberofvertices; + index=md.mesh.elements; + x=md.mesh.x; y=md.mesh.y; +else + numberofelements=md.mesh.numberofelements2d; + numberofnodes=md.mesh.numberofvertices2d; + index=md.mesh.elements2d; + x=md.mesh.x2d; y=md.mesh.y2d; +end + +%compute nodal functions coefficients N(x,y)=alpha x + beta y + gamma +[alpha beta]=GetNodalFunctionsCoeff(index,x,y); + +summation=[1;1;1]; +sx=(md.geometry.surface(index).*alpha)*summation; +sy=(md.geometry.surface(index).*beta)*summation; +s=sqrt(sx.^2+sy.^2); + +if md.mesh.dimension==3, + sx=project3d(md,'vector',sx,'type','element'); + sy=project3d(md,'vector',sy,'type','element'); + s=sqrt(sx.^2+sy.^2); +end Index: /issm/trunk-jpl/src/m/interp/ProfileValues.m =================================================================== --- /issm/trunk-jpl/src/m/interp/ProfileValues.m (revision 13008) +++ /issm/trunk-jpl/src/m/interp/ProfileValues.m (revision 13008) @@ -0,0 +1,21 @@ +function [Z,data_interp]=ProfileValues(md,data,xprof,yprof,resolution) +%PROFILEVALUES - compute the value of a field on a vertical profile +% +% This routine gets the value of a given field of the model on points +% given by filname (Argus type file) +% +% Usage: +% [z,data]=ProfileValues(md,data,filename,resolution) +% [z,data]=ProfileValues(md,data,profile_structure,resolution) + +%Get bed and surface for each 2d point, offset to make sure that it is inside the glacier system +offset=10^-3; +bed=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.bed,1),xprof,yprof)+offset; +surface=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.surface,1),xprof,yprof)-offset; + +%Some useful parameters +layers=ceil(mean(md.geometry.thickness)/resolution); +Z=(bed:resolution:surface)'; +X=xprof*ones(size(Z)); +Y=yprof*ones(size(Z)); +data_interp=InterpFromMeshToMesh3d(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.z,data,X,Y,Z,NaN); Index: /issm/trunk-jpl/src/m/interp/SectionValues.m =================================================================== --- /issm/trunk-jpl/src/m/interp/SectionValues.m (revision 13008) +++ /issm/trunk-jpl/src/m/interp/SectionValues.m (revision 13008) @@ -0,0 +1,127 @@ +function [index,X,Y,Z,S,data_interp]=SectionValues(md,data,infile,resolution) +%SECTIONVALUES - compute the value of a field on a section +% +% This routine gets the value of a given field of the model on points +% given by filname (Argus type file) +% +% Usage: +% [elements,x,y,z,s,data]=SectionValues(md,data,filename,resolution) +% [elements,x,y,z,s,data]=SectionValues(md,data,profile_structure,resolution) + +%check what we have for profile as input +if ischar(infile), + %read infile: + profile=expread(infile); + nods=profile.nods; + x=profile.x; + y=profile.y; +else + %read infile: + nods=infile.nods; + x=infile.x; + y=infile.y; +end + + +%get the specified resolution +if isnumeric(resolution(1)) + res_h=resolution(1); +else + error('SectionValues error message: wrong resolution type. Resolution must be an array [horizontal_resolution vertical_resolution]') +end +if md.mesh.dimension==3 + if (length(resolution)==2 & isnumeric(resolution(2))) + res_v=resolution(2); + else + error('SectionValues error message: wrong resolution type. Resolution must be an array [horizontal_resolution vertical_resolution]') + end +end + +%initialization +X=[]; %X-coordinate +Y=[]; %Y-coordinate +S=0; %curvilinear coordinate + +for i=1:nods-1 + + x_start=x(i); + x_end=x(i+1); + y_start=y(i); + y_end=y(i+1); + s_start=S(end); + + length_segment=sqrt((x_end-x_start)^2+(y_end-y_start)^2); + portion=ceil(length_segment/res_h); + + x_segment=zeros(portion,1); + y_segment=zeros(portion,1); + s_segment=zeros(portion,1); + + for j=1:portion + x_segment(j)=x_start+(j-1)*(x_end-x_start)/portion; + y_segment(j)=y_start+(j-1)*(y_end-y_start)/portion; + s_segment(j)=s_start+j*length_segment/portion; + end + + %plug into X and Y + X=[X;x_segment]; + Y=[Y;y_segment]; + S=[S;s_segment]; +end +X(end+1)=x(nods); +Y(end+1)=y(nods); + +%Number of nodes: +numberofnodes=size(X,1); + +%Compute Z +Z=zeros(numberofnodes,1); + +%New mesh and Data interpolation +if (md.mesh.dimension==2) + + %Interpolation of data on specified points + data_interp=InterpFromMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,data,X,Y); + %data_interp=InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,data,X,Y); + %data_interp=griddata(md.mesh.x,md.mesh.y,data,X,Y); + + %Compute index + index=[1:1:(numberofnodes-1);2:1:numberofnodes]'; + +else + + %vertically extrude mesh + + %Get bed and surface for each 2d point, offset to make sure that it is inside the glacier system + offset=10^-3; + bed=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.bed,1),X,Y)+offset; + surface=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.surface,1),X,Y)-offset; + + %Some useful parameters + layers=ceil(mean(md.geometry.thickness)/res_v); + nodesperlayer=numberofnodes; + nodestot=nodesperlayer*layers; + elementsperlayer=nodesperlayer-1; + elementstot=(nodesperlayer-1)*(layers-1); + + %initialization + X3=zeros(nodesperlayer*layers,1); Y3=zeros(nodesperlayer*layers,1); Z3=zeros(nodesperlayer*layers,1); S3=zeros(nodesperlayer*layers,1); index3=zeros(elementstot,4); + + %Get new coordinates in 3d + for i=1:layers + X3(i:layers:end)=X; + Y3(i:layers:end)=Y; + Z3(i:layers:end)=bed+(i-1)*(surface-bed)/(layers-1); + S3(i:layers:end)=S; + + if i=1)); - - -%average x and y over elements: -x_elem=md.mesh.x(md.mesh.elements)*[1;1;1]/3; -y_elem=md.mesh.y(md.mesh.elements)*[1;1;1]/3; - -while 1, - - %keep copy of new_flags for this loop: - new_flags_bak=new_flags; - - %extend new flags by connectivity - pos=find(new_flags); - - connected_elements=md.mesh.elementconnectivity(pos,:); - connected_elements=connected_elements(find(connected_elements)); - new_flags(connected_elements)=1; - - %get new elements: - new_elements=find(new_flags & ~new_flags_bak); - if ~length(new_elements), - %we are done! - break; - end - - - %check which of these new elements are more than distance away from the border elements - for i=1:length(new_elements), - dist=sqrt( (x_elem(border_elements)-x_elem(new_elements(i))).^2 + (y_elem(border_elements)-y_elem(new_elements(i))).^2)-distance; - if ~any(dist<0) - %none of the border elements are within distance, this element is outside out area of interest. - %ensure this element never gets found again in the connectivity. - pos=find(md.mesh.elementconnectivity==new_elements(i)); - md.mesh.elementconnectivity(pos)=0; - %exclude this new element from the new_flags! - new_flags(new_elements(i))=0; - end - end -end Index: sm/trunk-jpl/src/m/model/README =================================================================== --- /issm/trunk-jpl/src/m/model/README (revision 13007) +++ (revision ) @@ -1,11 +1,0 @@ -This directory is similar to @model, in that it only deals -with methods proper to the @model class. But we need a different -directory from @model, to make those methods public, not private. -The advantage of public methods is that they will use subsref and -susasgn to access data in a @model object. - -For ex: -calling md.x in a "@model" routine just accesses the x field in the md -structure. But calling md.x in a "model" routine will access the x field -through the susref routine in @model. This ensures that we protect the -data in @model classes from improper use. Index: sm/trunk-jpl/src/m/model/SectionValues.m =================================================================== --- /issm/trunk-jpl/src/m/model/SectionValues.m (revision 13007) +++ (revision ) @@ -1,127 +1,0 @@ -function [index,X,Y,Z,S,data_interp]=SectionValues(md,data,infile,resolution) -%SECTIONVALUES - compute the value of a field on a section -% -% This routine gets the value of a given field of the model on points -% given by filname (Argus type file) -% -% Usage: -% [elements,x,y,z,s,data]=SectionValues(md,data,filename,resolution) -% [elements,x,y,z,s,data]=SectionValues(md,data,profile_structure,resolution) - -%check what we have for profile as input -if ischar(infile), - %read infile: - profile=expread(infile); - nods=profile.nods; - x=profile.x; - y=profile.y; -else - %read infile: - nods=infile.nods; - x=infile.x; - y=infile.y; -end - - -%get the specified resolution -if isnumeric(resolution(1)) - res_h=resolution(1); -else - error('SectionValues error message: wrong resolution type. Resolution must be an array [horizontal_resolution vertical_resolution]') -end -if md.mesh.dimension==3 - if (length(resolution)==2 & isnumeric(resolution(2))) - res_v=resolution(2); - else - error('SectionValues error message: wrong resolution type. Resolution must be an array [horizontal_resolution vertical_resolution]') - end -end - -%initialization -X=[]; %X-coordinate -Y=[]; %Y-coordinate -S=0; %curvilinear coordinate - -for i=1:nods-1 - - x_start=x(i); - x_end=x(i+1); - y_start=y(i); - y_end=y(i+1); - s_start=S(end); - - length_segment=sqrt((x_end-x_start)^2+(y_end-y_start)^2); - portion=ceil(length_segment/res_h); - - x_segment=zeros(portion,1); - y_segment=zeros(portion,1); - s_segment=zeros(portion,1); - - for j=1:portion - x_segment(j)=x_start+(j-1)*(x_end-x_start)/portion; - y_segment(j)=y_start+(j-1)*(y_end-y_start)/portion; - s_segment(j)=s_start+j*length_segment/portion; - end - - %plug into X and Y - X=[X;x_segment]; - Y=[Y;y_segment]; - S=[S;s_segment]; -end -X(end+1)=x(nods); -Y(end+1)=y(nods); - -%Number of nodes: -numberofnodes=size(X,1); - -%Compute Z -Z=zeros(numberofnodes,1); - -%New mesh and Data interpolation -if (md.mesh.dimension==2) - - %Interpolation of data on specified points - data_interp=InterpFromMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,data,X,Y); - %data_interp=InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,data,X,Y); - %data_interp=griddata(md.mesh.x,md.mesh.y,data,X,Y); - - %Compute index - index=[1:1:(numberofnodes-1);2:1:numberofnodes]'; - -else - - %vertically extrude mesh - - %Get bed and surface for each 2d point, offset to make sure that it is inside the glacier system - offset=10^-3; - bed=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.bed,1),X,Y)+offset; - surface=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.surface,1),X,Y)-offset; - - %Some useful parameters - layers=ceil(mean(md.geometry.thickness)/res_v); - nodesperlayer=numberofnodes; - nodestot=nodesperlayer*layers; - elementsperlayer=nodesperlayer-1; - elementstot=(nodesperlayer-1)*(layers-1); - - %initialization - X3=zeros(nodesperlayer*layers,1); Y3=zeros(nodesperlayer*layers,1); Z3=zeros(nodesperlayer*layers,1); S3=zeros(nodesperlayer*layers,1); index3=zeros(elementstot,4); - - %Get new coordinates in 3d - for i=1:layers - X3(i:layers:end)=X; - Y3(i:layers:end)=Y; - Z3(i:layers:end)=bed+(i-1)*(surface-bed)/(layers-1); - S3(i:layers:end)=S; - - if i50000); - -%check that there is a GL -if isempty(pos_GL) - error('ThicknessCorrection error message: no grounding line has been detected. Check the model mask'); -end - -%get distance -if nargin==2, - distance=varargin{1}; -else - distance=10000; -end - -%modify thickness -if (debug), fprintf('%s',' correction progress: 0.00 %'); end -for i=1:length(pos_shelf) - - if (debug & mod(i,100)==0), - fprintf('\b\b\b\b\b\b\b%5.2f%s',i/length(pos_shelf)*100,' %'); - end - - %search the node on ice sheet the closest to i - [d posd]=min(sqrt((md.mesh.x(pos_shelf(i))-md.mesh.x(pos_GL)).^2+(md.mesh.y(pos_shelf(i))-md.mesh.y(pos_GL)).^2)); - - if d>distance, - - %if d > 15km, hydrostatic equilibrium - hydrostatic_ratio(pos_shelf(i))=1; - thickness(pos_shelf(i))=thickness_hydro(pos_shelf(i)); - - else - - %else: quadratic combination of hydrostatic equilibrium and observations - hydrostatic_ratio(pos_shelf(i))=(d/distance)^2; - thickness(pos_shelf(i))=(1-hydrostatic_ratio(pos_shelf(i)))*thickness(pos_shelf(i))+hydrostatic_ratio(pos_shelf(i))*thickness_hydro(pos_shelf(i)); - - end -end -if (debug), fprintf('\b\b\b\b\b\b\b%5.2f%s\n',100,' %'); end - -%check the computed thickness -minth=1/(1-md.materials.rho_ice/md.materials.rho_water); -pos=find(isnan(thickness) | (thickness<=0)); -thickness(pos)=minth; -hydrostatic_ratio(pos)=-1; - -%change bed to take into account the changes in thickness -md.geometry.thickness=thickness; -md.geometry.hydrostatic_ratio=hydrostatic_ratio; -md.geometry.bed=md.geometry.surface-md.geometry.thickness; Index: /issm/trunk-jpl/src/m/model/inversions/misfit.m =================================================================== --- /issm/trunk-jpl/src/m/model/inversions/misfit.m (revision 13008) +++ /issm/trunk-jpl/src/m/model/inversions/misfit.m (revision 13008) @@ -0,0 +1,37 @@ +function J=misfit(md) +%MISFIT - compute misfit +% +% Usage: +% J=misfit(md) +% +% Example: +% J=misfit(md) +% + +if md.mesh.dimension==2, + elements=md.mesh.elements; + x=md.mesh.x; + y=md.mesh.y; + vx=md.initialization.vx; + vy=md.initialization.vy; + vx_obs=md.inversion.vx_obs; + vy_obs=md.inversion.vy_obs; +else + elements=md.mesh.elements2d; + x=md.mesh.x2d; + y=md.mesh.y2d; + vx=project2d(md,md.initialization.vx,md.mesh.numberoflayers); + vy=project2d(md,md.initialization.vy,md.mesh.numberoflayers); + vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers); + vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers); +end + +%compute areas; +areas=GetAreas(elements,x,y); + +%compute delta v on elements +deltav=1/2*( (vx-vx_obs).^2+(vy-vy_obs).^2)/md.constants.yts^2; +deltav_elem=deltav(elements)*[1;1;1]/3; + +%compute misfit +J=sum(deltav_elem.*areas); Index: /issm/trunk-jpl/src/m/model/mesh/bamg.m =================================================================== --- /issm/trunk-jpl/src/m/model/mesh/bamg.m (revision 13007) +++ /issm/trunk-jpl/src/m/model/mesh/bamg.m (revision 13008) @@ -339,2 +339,145 @@ error('Output mesh has orphans. Decrease MaxCornerAngle to prevent outside points (ex: 0.01)'); end +end + +function geom=processgeometry(geom,tol,outline); % {{{ + +%Deal with edges +disp('Checking Edge crossing...'); +i=0; +while (ii); + geom.Edges(posedges)=geom.Edges(posedges)-1; + + %update counter + i=i-1; + end +end +if num, + disp(['WARNING: ' num2str(num) ' points outside the domain outline have been removed']); +end + +%Check point spacing +if ~isnan(tol), + disp('Checking point spacing...'); + i=0; + while (ij); + geom.Edges(posedges)=geom.Edges(posedges)-1; + + %update counter + j=j-1; + + end + end + end +end +%remove empty edges +geom.Edges(find(geom.Edges(:,1)==geom.Edges(:,2)),:)=[]; +end % }}} Index: sm/trunk-jpl/src/m/model/misfit.m =================================================================== --- /issm/trunk-jpl/src/m/model/misfit.m (revision 13007) +++ (revision ) @@ -1,37 +1,0 @@ -function J=misfit(md) -%MISFIT - compute misfit -% -% Usage: -% J=misfit(md) -% -% Example: -% J=misfit(md) -% - -if md.mesh.dimension==2, - elements=md.mesh.elements; - x=md.mesh.x; - y=md.mesh.y; - vx=md.initialization.vx; - vy=md.initialization.vy; - vx_obs=md.inversion.vx_obs; - vy_obs=md.inversion.vy_obs; -else - elements=md.mesh.elements2d; - x=md.mesh.x2d; - y=md.mesh.y2d; - vx=project2d(md,md.initialization.vx,md.mesh.numberoflayers); - vy=project2d(md,md.initialization.vy,md.mesh.numberoflayers); - vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers); - vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers); -end - -%compute areas; -areas=GetAreas(elements,x,y); - -%compute delta v on elements -deltav=1/2*( (vx-vx_obs).^2+(vy-vy_obs).^2)/md.constants.yts^2; -deltav_elem=deltav(elements)*[1;1;1]/3; - -%compute misfit -J=sum(deltav_elem.*areas); Index: sm/trunk-jpl/src/m/model/modelsextract.m =================================================================== --- /issm/trunk-jpl/src/m/model/modelsextract.m (revision 13007) +++ (revision ) @@ -1,87 +1,0 @@ -function md_list=modelsextract(md,flags,minel,varargin) -%modelsextract - extract several self contained models according to a list of element flags. -% -% The difference between this routine and the modelextract.m routine (without an 's') is that -% as many models are extracted as there are closed contours defined in area. -% This routine is needed for example when doing data assimilation of ice shelves in Antarctica. -% Many independent ice shelves are present, and we don't want data assimilation on one ice shelf -% to be hindered by another totally independent ice shelf. -% -% Usage: -% md_list=modelsextract(md,elementfalgs,minel); -% -% Examples: -% md_list=modelsextract(md,md.mask.elementonfloatingice,1000); -% -% See also: EXTRUDE, COLLAPSE, MODELEXTRACT - -disp('selecting pools of elements'); -%go through flags and build as many independent element flags as there are groups of connected 1s -%in flags. - -%2D or 3D? -if md.mesh.dimension==3, - numberofelements=md.mesh.numberofelements2d; %this will be forgotten when we get out. - flags=project2d(md,flags,1); -else - numberofelements=md.mesh.numberofelements; -end - -%recover extra arguments: -distance=0; -if nargin==4, - distance=varargin{1}; -end - -flag_list=cell(0,1); - -for i=1:size(flags,1), - - if (flags(i)), - - %ok, we are sure element i is part of a new pool. - pool=zeros(numberofelements,1); - pool=PropagateFlagsFromConnectivity(md.mesh.elementconnectivity,pool,i,flags); - flag_list{end+1,1}=pool; - - %speed up rest of computation by taking pool out of flags: - pos=find(pool);flags(pos)=0; - - end -end - -%go through flag_list and discard any pool of less than minel elements: -ex_pos=[]; -for i=1:length(flag_list), - if length(find(flag_list{i}))0)))=1; Index: /issm/trunk-jpl/src/m/model/plot/radarpower.m =================================================================== --- /issm/trunk-jpl/src/m/model/plot/radarpower.m (revision 13008) +++ /issm/trunk-jpl/src/m/model/plot/radarpower.m (revision 13008) @@ -0,0 +1,118 @@ +function md=radarpower(md,varargin) +%RADARPOWER - overlay a power radar image on an existing mesh +% +% This routine will overlay a power radar image on an existing mesh. +% The power amplitude will be output to vel for now. +% In the future, think about a field to hold this value. +% +% Usage: +% md=radarpower(md,options); +% md=radarpower(md) + +%If gdal does not work, uncomment the following line +%setenv('LD_LIBRARY_PATH','/proj/ice/larour/issm/trunk/externalpackages/gdal/install/lib/'); +%Parse inputs +if nargin==1, + options=pairoptions; +else + options=varargin{:}; + if ~isa(options,'pairoptions'), + options=pairoptions(varargin{:}); + end +end + +highres=getfieldvalue(options,'highres',0); +xlim=getfieldvalue(options,'xlim',[min(md.mesh.x) max(md.mesh.x)]); +ylim=getfieldvalue(options,'ylim',[min(md.mesh.y) max(md.mesh.y)]); +posting=getfieldvalue(options,'posting',0); % 0 -> image posting default + +%find gdal coordinates +x0=min(xlim); x1=max(xlim); +y0=min(ylim); y1=max(ylim); + +%figure out if we should go look for Greenland or Antarctica geotiff, or if user provided one. +if ~exist(options,'overlay_image'), + if strcmpi(md.mesh.hemisphere,'n'), + if ~exist([jplsvn() '/projects/ModelData/MOG/mog150_greenland_map.jpg']), + error(['radarpower error message: file ' jplsvn() '/projects/ModelData/MOG/mog150_greenland_map.jpg not found.']); + end + name = 'mog150_greenland_map'; + %name = 'mog100_hp1_v10'; + %name = 'mog500_hp1_v10'; + jpgim=[jplsvn() '/projects/ModelData/MOG/' name '.jpg']; + geom=load([jplsvn() '/projects/ModelData/MOG/' name '.jpgw'],'ascii'); + + %geom: xposting nbcols nbrows yposting xmin ymax + xmin=max(geom(5),x0); + xmax=min(geom(5)+geom(1)*geom(2),x1); + ymin=max(geom(6)-geom(3)*geom(4),y0); + ymax=min(geom(6),y1); + + firstcol=max(1,floor((xmin-geom(5))/geom(1))); %x min + firstrow=max(1,floor((geom(6)-ymax)/geom(4))); %y max + numcols=floor((xmax-xmin)/geom(1)); % x posting + numrows=floor((ymax-ymin)/geom(4)); % y posting + pixelskip=max(1,ceil(posting/geom(1))); + + %Read and crop file + disp('Warning: expecting coordinates in polar stereographic (Std Latitude: 70ºN Meridian: 45º)'); + im=imread(jpgim); + im=im(firstrow:firstrow+numrows-1,firstcol:firstcol+numcols-1); + md.radaroverlay.pwr=double(flipud(im(1:pixelskip:end,1:pixelskip:end))); + md.radaroverlay.x=(xmin:(xmax-xmin)/(size(md.radaroverlay.pwr,2)-1):xmax); + md.radaroverlay.y=(ymin:(ymax-ymin)/(size(md.radaroverlay.pwr,1)-1):ymax); + + elseif strcmpi(md.mesh.hemisphere,'s'), + if highres, + if ~exist([jplsvn() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif']), + error(['radarpower error message: file ' jplsvn() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif not found.']); + end + geotiff_name=[jplsvn() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif']; + else + if ~exist([jplsvn() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif']), + error(['radarpower error message: file ' jplsvn() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif not found.']); + end + geotiff_name=[jplsvn() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif']; + end + + %Name of image + inputname='./temp.tif'; + eval(['!gdal_translate -quiet -projwin ' num2str(x0) ' ' num2str(y1) ' ' num2str(x1) ' ' num2str(y0) ' ' geotiff_name ' ' inputname ]); + + %Read in temp.tif: + im=imread('temp.tif','TIFF'); + pixelskip=max(1,ceil(posting/((x1-x0)/(size(im,2))))); + md.radaroverlay.pwr=double(flipud(im(1:pixelskip:end,1:pixelskip:end))); + md.radaroverlay.x=(x0:(x1-x0)/(size(md.radaroverlay.pwr,2)-1):x1); + md.radaroverlay.y=(y0:(y1-y0)/(size(md.radaroverlay.pwr,1)-1):y1); + + %Erase image + system('rm -rf ./temp.tif'); + + else + error('field hemisphere should either be ''n'' or ''s'''); + end +else + %ok, user provided an image. check we also have overlay_xlim and overlay_ylim options, to know what range of coordinates the image covers. + if (~exist(options,'overlay_xlim') | ~exist(options,'overlay_xlim')| ~exist(options,'overlay_xposting')| ~exist(options,'overlay_yposting')), + error('radarpower error message: please provide overlay_xlim, overlay_ylim, overlay_xposting and overlay_yposting options together with overlay_image option'); + end + overlay_image=getfieldvalue(options,'overlay_image'); + overlay_xlim=getfieldvalue(options,'overlay_xlim'); + overlay_ylim=getfieldvalue(options,'overlay_ylim'); + overlay_xposting=getfieldvalue(options,'overlay_xposting'); + overlay_yposting=getfieldvalue(options,'overlay_yposting'); + + sizex=floor((x1-x0)/overlay_xposting); + sizey=floor((y1-y0)/overlay_yposting); + topleftx=floor((x0-overlay_xlim(1))/overlay_xposting); % x min + toplefty=floor((overlay_ylim(2)-y1)/overlay_yposting); % y max + + %Read and crop file + disp('Warning: expecting coordinates in polar stereographic (Std Latitude: 70ºN Meridian: 45º)'); + im=imread(overlay_image); + im=im(toplefty:toplefty+sizey,topleftx:topleftx+sizex); + md.radaroverlay.pwr=double(flipud(im)); + md.radaroverlay.x=(x0:(x1-x0)/(size(md.radaroverlay.pwr,2)-1):x1); + md.radaroverlay.y=(y0:(y1-y0)/(size(md.radaroverlay.pwr,1)-1):y1); +end Index: sm/trunk-jpl/src/m/model/processgeometry.m =================================================================== --- /issm/trunk-jpl/src/m/model/processgeometry.m (revision 13007) +++ (revision ) @@ -1,140 +1,0 @@ -function geom=processgeometry(geom,tol,outline); - -%Deal with edges -disp('Checking Edge crossing...'); -i=0; -while (ii); - geom.Edges(posedges)=geom.Edges(posedges)-1; - - %update counter - i=i-1; - end -end -if num, - disp(['WARNING: ' num2str(num) ' points outside the domain outline have been removed']); -end - -%Check point spacing -if ~isnan(tol), - disp('Checking point spacing...'); - i=0; - while (ij); - geom.Edges(posedges)=geom.Edges(posedges)-1; - - %update counter - j=j-1; - - end - end - end -end -%remove empty edges -geom.Edges(find(geom.Edges(:,1)==geom.Edges(:,2)),:)=[]; Index: sm/trunk-jpl/src/m/model/qstat.m =================================================================== --- /issm/trunk-jpl/src/m/model/qstat.m (revision 13007) +++ (revision ) @@ -1,8 +1,0 @@ -function qstat(md) -%QSTAT - check job status on remote cluster -% -% Usage: -% qstat(md) - -%run qstat command on remote cluster -issmssh(md.cluster,['qstat -a']); Index: sm/trunk-jpl/src/m/model/radarpower.m =================================================================== --- /issm/trunk-jpl/src/m/model/radarpower.m (revision 13007) +++ (revision ) @@ -1,118 +1,0 @@ -function md=radarpower(md,varargin) -%RADARPOWER - overlay a power radar image on an existing mesh -% -% This routine will overlay a power radar image on an existing mesh. -% The power amplitude will be output to vel for now. -% In the future, think about a field to hold this value. -% -% Usage: -% md=radarpower(md,options); -% md=radarpower(md) - -%If gdal does not work, uncomment the following line -%setenv('LD_LIBRARY_PATH','/proj/ice/larour/issm/trunk/externalpackages/gdal/install/lib/'); -%Parse inputs -if nargin==1, - options=pairoptions; -else - options=varargin{:}; - if ~isa(options,'pairoptions'), - options=pairoptions(varargin{:}); - end -end - -highres=getfieldvalue(options,'highres',0); -xlim=getfieldvalue(options,'xlim',[min(md.mesh.x) max(md.mesh.x)]); -ylim=getfieldvalue(options,'ylim',[min(md.mesh.y) max(md.mesh.y)]); -posting=getfieldvalue(options,'posting',0); % 0 -> image posting default - -%find gdal coordinates -x0=min(xlim); x1=max(xlim); -y0=min(ylim); y1=max(ylim); - -%figure out if we should go look for Greenland or Antarctica geotiff, or if user provided one. -if ~exist(options,'overlay_image'), - if strcmpi(md.mesh.hemisphere,'n'), - if ~exist([jplsvn() '/projects/ModelData/MOG/mog150_greenland_map.jpg']), - error(['radarpower error message: file ' jplsvn() '/projects/ModelData/MOG/mog150_greenland_map.jpg not found.']); - end - name = 'mog150_greenland_map'; - %name = 'mog100_hp1_v10'; - %name = 'mog500_hp1_v10'; - jpgim=[jplsvn() '/projects/ModelData/MOG/' name '.jpg']; - geom=load([jplsvn() '/projects/ModelData/MOG/' name '.jpgw'],'ascii'); - - %geom: xposting nbcols nbrows yposting xmin ymax - xmin=max(geom(5),x0); - xmax=min(geom(5)+geom(1)*geom(2),x1); - ymin=max(geom(6)-geom(3)*geom(4),y0); - ymax=min(geom(6),y1); - - firstcol=max(1,floor((xmin-geom(5))/geom(1))); %x min - firstrow=max(1,floor((geom(6)-ymax)/geom(4))); %y max - numcols=floor((xmax-xmin)/geom(1)); % x posting - numrows=floor((ymax-ymin)/geom(4)); % y posting - pixelskip=max(1,ceil(posting/geom(1))); - - %Read and crop file - disp('Warning: expecting coordinates in polar stereographic (Std Latitude: 70ºN Meridian: 45º)'); - im=imread(jpgim); - im=im(firstrow:firstrow+numrows-1,firstcol:firstcol+numcols-1); - md.radaroverlay.pwr=double(flipud(im(1:pixelskip:end,1:pixelskip:end))); - md.radaroverlay.x=(xmin:(xmax-xmin)/(size(md.radaroverlay.pwr,2)-1):xmax); - md.radaroverlay.y=(ymin:(ymax-ymin)/(size(md.radaroverlay.pwr,1)-1):ymax); - - elseif strcmpi(md.mesh.hemisphere,'s'), - if highres, - if ~exist([jplsvn() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif']), - error(['radarpower error message: file ' jplsvn() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif not found.']); - end - geotiff_name=[jplsvn() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif']; - else - if ~exist([jplsvn() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif']), - error(['radarpower error message: file ' jplsvn() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif not found.']); - end - geotiff_name=[jplsvn() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif']; - end - - %Name of image - inputname='./temp.tif'; - eval(['!gdal_translate -quiet -projwin ' num2str(x0) ' ' num2str(y1) ' ' num2str(x1) ' ' num2str(y0) ' ' geotiff_name ' ' inputname ]); - - %Read in temp.tif: - im=imread('temp.tif','TIFF'); - pixelskip=max(1,ceil(posting/((x1-x0)/(size(im,2))))); - md.radaroverlay.pwr=double(flipud(im(1:pixelskip:end,1:pixelskip:end))); - md.radaroverlay.x=(x0:(x1-x0)/(size(md.radaroverlay.pwr,2)-1):x1); - md.radaroverlay.y=(y0:(y1-y0)/(size(md.radaroverlay.pwr,1)-1):y1); - - %Erase image - system('rm -rf ./temp.tif'); - - else - error('field hemisphere should either be ''n'' or ''s'''); - end -else - %ok, user provided an image. check we also have overlay_xlim and overlay_ylim options, to know what range of coordinates the image covers. - if (~exist(options,'overlay_xlim') | ~exist(options,'overlay_xlim')| ~exist(options,'overlay_xposting')| ~exist(options,'overlay_yposting')), - error('radarpower error message: please provide overlay_xlim, overlay_ylim, overlay_xposting and overlay_yposting options together with overlay_image option'); - end - overlay_image=getfieldvalue(options,'overlay_image'); - overlay_xlim=getfieldvalue(options,'overlay_xlim'); - overlay_ylim=getfieldvalue(options,'overlay_ylim'); - overlay_xposting=getfieldvalue(options,'overlay_xposting'); - overlay_yposting=getfieldvalue(options,'overlay_yposting'); - - sizex=floor((x1-x0)/overlay_xposting); - sizey=floor((y1-y0)/overlay_yposting); - topleftx=floor((x0-overlay_xlim(1))/overlay_xposting); % x min - toplefty=floor((overlay_ylim(2)-y1)/overlay_yposting); % y max - - %Read and crop file - disp('Warning: expecting coordinates in polar stereographic (Std Latitude: 70ºN Meridian: 45º)'); - im=imread(overlay_image); - im=im(toplefty:toplefty+sizey,topleftx:topleftx+sizex); - md.radaroverlay.pwr=double(flipud(im)); - md.radaroverlay.x=(x0:(x1-x0)/(size(md.radaroverlay.pwr,2)-1):x1); - md.radaroverlay.y=(y0:(y1-y0)/(size(md.radaroverlay.pwr,1)-1):y1); -end Index: sm/trunk-jpl/src/m/model/recover_areas.m =================================================================== --- /issm/trunk-jpl/src/m/model/recover_areas.m (revision 13007) +++ (revision ) @@ -1,22 +1,0 @@ -function [hutterflag macayealflag pattynflag stokesflag filltype]=recover_areas(md,varargin); -%RECOVER_AREAS - flag the element depending on the physical model that is assigned to them -% -% This routine is called by setelementstype, do not use -% -% Usage: -% [hutterflag macayealflag pattynflag stokesflag filltype]=recover_areas(md,varargin); - - %go through varargin, extract options and plug them into subtype options, by order of appearance - options=pairoptions(varargin{:}); - options=deleteduplicates(options,1); - - %recover elements distribution - hutterflag =FlagElements(md,getfieldvalue(options,'hutter','')); - macayealflag=FlagElements(md,getfieldvalue(options,'macayeal','')); - pattynflag =FlagElements(md,getfieldvalue(options,'pattyn','')); - stokesflag =FlagElements(md,getfieldvalue(options,'stokes','')); - filltype =getfieldvalue(options,'fill','none'); - -end %end function - - Index: sm/trunk-jpl/src/m/model/recover_areas.py =================================================================== --- /issm/trunk-jpl/src/m/model/recover_areas.py (revision 13007) +++ (revision ) @@ -1,26 +1,0 @@ -from pairoptions import * -from FlagElements import * - -def recover_areas(md,*args): - """ - RECOVER_AREAS - flag the element depending on the physical model that is assigned to them - - This routine is called by setelementstype, do not use - - Usage: - [hutterflag macayealflag pattynflag stokesflag filltype]=recover_areas(md,varargin); - """ - - #go through varargin, extract options and plug them into subtype options, by order of appearance - options=pairoptions(*args) -# options=deleteduplicates(options,1); - - #recover elements distribution - hutterflag =FlagElements(md,options.getfieldvalue('hutter','')) - macayealflag=FlagElements(md,options.getfieldvalue('macayeal','')) - pattynflag =FlagElements(md,options.getfieldvalue('pattyn','')) - stokesflag =FlagElements(md,options.getfieldvalue('stokes','')) - filltype =options.getfieldvalue('fill','none') - - return hutterflag,macayealflag,pattynflag,stokesflag,filltype - Index: /issm/trunk-jpl/src/m/model/regional/regionaltransient2d.m =================================================================== --- /issm/trunk-jpl/src/m/model/regional/regionaltransient2d.m (revision 13008) +++ /issm/trunk-jpl/src/m/model/regional/regionaltransient2d.m (revision 13008) @@ -0,0 +1,158 @@ +function md2=regionaltransient2d(md1,area,hmin,hmax,err,stepres) +%regionaltransient2d - extract a model according to an Argus contour or flag list and remesh +% at new resolution res +% +% This routine extracts a submodel from a bigger model with respect to a given contour +% md must be followed by the corresponding exp domain file (argus type, .exp extension). +% The model will be remeshed at high rsolution hmin and low resolution hmax. The ice +% boundary velocities will be spc'd to the transient velocities at saved transient steps +% at the resolution optionally provided for stepres. A stepres of 2 means that you wish +% to skip every other saved transient step. This is useful when extracting a long transient. +% +% Usage: +% md2=regionaltransient2d(md1,area,hmin,hmax,err); +% +% Examples: +% md2=regionaltransient2d(md,'Domain.exp',500,10000,[15 250]); +% md2=regionaltransient2d(md,'Domain.exp',3000,15000,[10 300],2); +% +% See also: MODELEXTRACT, EXTRUDE, COLLAPSE + +%some checks +if ((nargin~=5) & (nargin~=6)), + help regionaltransient2d + error('regionaltransient2d error message: bad usage'); +end + +%get check option +if (nargin==5), + stepres=1; +end + +%take every fields from model +mde=modelextract(md1,area); +mde.private.bamg=[]; +mde.mesh.extractedvertices=nan; +mde.mesh.extractedelements=nan; + +%remesh +md2=bamg(mde,'hmin',hmin,'hmax',hmax,'field',[mde.inversion.vel_obs mde.geometry.surface],'splitcorner',1,'KeepVertices',0,'err',err); +md2=setmask(md2,'',''); + +%automatically modify fields + + %loop over model fields + model_fields=fields(md1); + for i=1:length(model_fields), + + %get field + field=md1.(model_fields{i}); + fieldsize=size(field); + + %copy field, interpolated to new mesh + if isobject(field), %recursive call + object_fields=fields(md1.(model_fields{i})); + fname=['(model_fields{i}).(object_fields{j})']; + else + object_fields=field; + fname=['(model_fields{i})']; + end + for j=1:length(object_fields), + %get field + field=eval(['md2.' fname]); + fieldsize=size(field); + + %size = number of nodes * n + for n=1:fieldsize(2) + if fieldsize(1)==mde.mesh.numberofvertices + if(sum(field(:,n) ~= field(1,n)) == 0) + eval(['md2.' fname '(1:md2.mesh.numberofvertices,n)=field(1,n)*ones(md2.mesh.numberofvertices,1);']); + else + eval(['md2.' fname '(1:md2.mesh.numberofvertices,n)=InterpFromMeshToMesh2d(mde.mesh.elements,mde.mesh.x,mde.mesh.y,field(:,n),md2.mesh.x,md2.mesh.y);']); + end + eval(['md2.' fname '(:,n)=md2.' fname '(1:md2.mesh.numberofvertices,n);']); + elseif fieldsize(1)==mde.mesh.numberofvertices+1 + if(sum(field(1:end-1,n) ~= field(1,n)) == 0) + eval(['md2.' fname '(1:md2.mesh.numberofvertices+1,n)=[field(1,n)*ones(md2.mesh.numberofvertices,1); field(end,n)];']); + else + eval(['md2.' fname '(1:md2.mesh.numberofvertices+1,n)=[InterpFromMeshToMesh2d(mde.mesh.elements,mde.mesh.x,mde.mesh.y,field(1:end-1,n),md2.mesh.x,md2.mesh.y); field(end,n)];']); + end + eval(['md2.' fname '(:,n)=md2.' fname '(1:md2.mesh.numberofvertices+1,n)']); + %size = number of elements * n + elseif fieldsize(1)==mde.mesh.numberofelements + if(sum(field(1:end-1,n) ~= field(1,n)) == 0) + eval(['md2.' fname '(1:md2.mesh.numberofelements,n)=field(1,n)*ones(md2.mesh.numberofelements,1);']); + else + eval(['md2.' fname '(1:md2.mesh.numberofelements,n)=InterpFromMeshToMesh2d(mde.mesh.elements,mde.mesh.x,mde.mesh.y,field(:,n),md2.mesh.x,md2.mesh.y);']); + end + eval(['md2.' fname '(:,n)=md2.' fname '(1:md2.mesh.numberofelements,n);']); + end + end + end + end + + %Read transient velocities and thickness, looping through only the populated times + spcx=[]; + spcy=[]; + spct=[]; + steps=[]; + nsteps=length(md1.results.TransientSolution); + count=0; + numElements=arrayfun(@(x) numel(x.step), md1.results.TransientSolution); + for t=find(numElements==1) + if ~isempty(md1.results.TransientSolution(t).Vel) & mod(count,stepres)==0, + vx=PatchToVec(md1.results.TransientSolution(t).Vx); + vy=PatchToVec(md1.results.TransientSolution(t).Vy); + thickness=PatchToVec(md1.results.TransientSolution(t).Thickness); + spcx=[spcx InterpFromMeshToMesh2d(md1.mesh.elements,md1.mesh.x,md1.mesh.y,vx,md2.mesh.x,md2.mesh.y)]; + spcy=[spcy InterpFromMeshToMesh2d(md1.mesh.elements,md1.mesh.x,md1.mesh.y,vy,md2.mesh.x,md2.mesh.y)]; + spct=[spct InterpFromMeshToMesh2d(md1.mesh.elements,md1.mesh.x,md1.mesh.y,thickness,md2.mesh.x,md2.mesh.y)]; + steps=[steps t*md1.timestepping.time_step]; + end + count=count+1; + end + + %As long as there are recorded time steps, spc the boundaries with velocities + if nsteps > 0 + md2.diagnostic.spcvx=md2.diagnostic.spcvx*ones(1,size(spcx,2)); + md2.diagnostic.spcvy=md2.diagnostic.spcvy*ones(1,size(spcy,2)); + md2.diagnostic.spcvz=md2.diagnostic.spcvz*ones(1,size(spcx,2)); + md2.prognostic.spcthickness=md2.prognostic.spcthickness*ones(1,size(spct,2)); + md2.diagnostic.spcvx(find(md2.mesh.vertexonboundary),:)=spcx(find(md2.mesh.vertexonboundary),:); + md2.diagnostic.spcvy(find(md2.mesh.vertexonboundary),:)=spcy(find(md2.mesh.vertexonboundary),:); + md2.diagnostic.spcvz(find(md2.mesh.vertexonboundary),:)=0; + md2.prognostic.spcthickness(find(md2.mesh.vertexonboundary),:)=spct(find(md2.mesh.vertexonboundary),:); + md2.diagnostic.spcvx=[md2.diagnostic.spcvx; steps]; + md2.diagnostic.spcvy=[md2.diagnostic.spcvy; steps]; + md2.diagnostic.spcvz=[md2.diagnostic.spcvz; steps]; + md2.prognostic.spcthickness=[md2.prognostic.spcthickness; steps]; + end + + %Diagnostic. Don't spc the icefront vertices. + if ~isnan(md2.diagnostic.icefront) + md1s=modelextract(md1,area); + %md2.diagnostic.icefront=[md2.mesh.segments 2]; + e2=md2.mesh.segments(:,end); + e1=md1s.mesh.segments(:,end); + + pload = nan*ones(size(md1s.mesh.elements,1),1); + pload(md1s.diagnostic.icefront(:,end-1))=md1s.diagnostic.icefront(:,end); + + x2=mean(md2.mesh.x(md2.mesh.elements(e2,:)),2); + y2=mean(md2.mesh.y(md2.mesh.elements(e2,:)),2); + x1=mean(md1s.mesh.x(md1s.mesh.elements),2); + y1=mean(md1s.mesh.y(md1s.mesh.elements),2); + + pload2=griddata(x1,y1,pload,x2,y2,'nearest'); + md2.diagnostic.icefront=[md2.mesh.segments(~isnan(pload2),:) pload2(~isnan(pload2))]; + md2.diagnostic.spcvx(unique(md2.diagnostic.icefront(:,1:2)),:)=nan; + md2.diagnostic.spcvy(unique(md2.diagnostic.icefront(:,1:2)),:)=nan; + md2.diagnostic.spcvz(unique(md2.diagnostic.icefront(:,1:2)),:)=nan; + md2.prognostic.spcthickness(unique(md2.diagnostic.icefront(:,1:2)),:)=nan; + end + + %Clear results fields + if isstruct(md1.results), + md2.results=[]; + end + Index: sm/trunk-jpl/src/m/model/regionaltransient2d.m =================================================================== --- /issm/trunk-jpl/src/m/model/regionaltransient2d.m (revision 13007) +++ (revision ) @@ -1,158 +1,0 @@ -function md2=regionaltransient2d(md1,area,hmin,hmax,err,stepres) -%regionaltransient2d - extract a model according to an Argus contour or flag list and remesh -% at new resolution res -% -% This routine extracts a submodel from a bigger model with respect to a given contour -% md must be followed by the corresponding exp domain file (argus type, .exp extension). -% The model will be remeshed at high rsolution hmin and low resolution hmax. The ice -% boundary velocities will be spc'd to the transient velocities at saved transient steps -% at the resolution optionally provided for stepres. A stepres of 2 means that you wish -% to skip every other saved transient step. This is useful when extracting a long transient. -% -% Usage: -% md2=regionaltransient2d(md1,area,hmin,hmax,err); -% -% Examples: -% md2=regionaltransient2d(md,'Domain.exp',500,10000,[15 250]); -% md2=regionaltransient2d(md,'Domain.exp',3000,15000,[10 300],2); -% -% See also: MODELEXTRACT, EXTRUDE, COLLAPSE - -%some checks -if ((nargin~=5) & (nargin~=6)), - help regionaltransient2d - error('regionaltransient2d error message: bad usage'); -end - -%get check option -if (nargin==5), - stepres=1; -end - -%take every fields from model -mde=modelextract(md1,area); -mde.private.bamg=[]; -mde.mesh.extractedvertices=nan; -mde.mesh.extractedelements=nan; - -%remesh -md2=bamg(mde,'hmin',hmin,'hmax',hmax,'field',[mde.inversion.vel_obs mde.geometry.surface],'splitcorner',1,'KeepVertices',0,'err',err); -md2=setmask(md2,'',''); - -%automatically modify fields - - %loop over model fields - model_fields=fields(md1); - for i=1:length(model_fields), - - %get field - field=md1.(model_fields{i}); - fieldsize=size(field); - - %copy field, interpolated to new mesh - if isobject(field), %recursive call - object_fields=fields(md1.(model_fields{i})); - fname=['(model_fields{i}).(object_fields{j})']; - else - object_fields=field; - fname=['(model_fields{i})']; - end - for j=1:length(object_fields), - %get field - field=eval(['md2.' fname]); - fieldsize=size(field); - - %size = number of nodes * n - for n=1:fieldsize(2) - if fieldsize(1)==mde.mesh.numberofvertices - if(sum(field(:,n) ~= field(1,n)) == 0) - eval(['md2.' fname '(1:md2.mesh.numberofvertices,n)=field(1,n)*ones(md2.mesh.numberofvertices,1);']); - else - eval(['md2.' fname '(1:md2.mesh.numberofvertices,n)=InterpFromMeshToMesh2d(mde.mesh.elements,mde.mesh.x,mde.mesh.y,field(:,n),md2.mesh.x,md2.mesh.y);']); - end - eval(['md2.' fname '(:,n)=md2.' fname '(1:md2.mesh.numberofvertices,n);']); - elseif fieldsize(1)==mde.mesh.numberofvertices+1 - if(sum(field(1:end-1,n) ~= field(1,n)) == 0) - eval(['md2.' fname '(1:md2.mesh.numberofvertices+1,n)=[field(1,n)*ones(md2.mesh.numberofvertices,1); field(end,n)];']); - else - eval(['md2.' fname '(1:md2.mesh.numberofvertices+1,n)=[InterpFromMeshToMesh2d(mde.mesh.elements,mde.mesh.x,mde.mesh.y,field(1:end-1,n),md2.mesh.x,md2.mesh.y); field(end,n)];']); - end - eval(['md2.' fname '(:,n)=md2.' fname '(1:md2.mesh.numberofvertices+1,n)']); - %size = number of elements * n - elseif fieldsize(1)==mde.mesh.numberofelements - if(sum(field(1:end-1,n) ~= field(1,n)) == 0) - eval(['md2.' fname '(1:md2.mesh.numberofelements,n)=field(1,n)*ones(md2.mesh.numberofelements,1);']); - else - eval(['md2.' fname '(1:md2.mesh.numberofelements,n)=InterpFromMeshToMesh2d(mde.mesh.elements,mde.mesh.x,mde.mesh.y,field(:,n),md2.mesh.x,md2.mesh.y);']); - end - eval(['md2.' fname '(:,n)=md2.' fname '(1:md2.mesh.numberofelements,n);']); - end - end - end - end - - %Read transient velocities and thickness, looping through only the populated times - spcx=[]; - spcy=[]; - spct=[]; - steps=[]; - nsteps=length(md1.results.TransientSolution); - count=0; - numElements=arrayfun(@(x) numel(x.step), md1.results.TransientSolution); - for t=find(numElements==1) - if ~isempty(md1.results.TransientSolution(t).Vel) & mod(count,stepres)==0, - vx=PatchToVec(md1.results.TransientSolution(t).Vx); - vy=PatchToVec(md1.results.TransientSolution(t).Vy); - thickness=PatchToVec(md1.results.TransientSolution(t).Thickness); - spcx=[spcx InterpFromMeshToMesh2d(md1.mesh.elements,md1.mesh.x,md1.mesh.y,vx,md2.mesh.x,md2.mesh.y)]; - spcy=[spcy InterpFromMeshToMesh2d(md1.mesh.elements,md1.mesh.x,md1.mesh.y,vy,md2.mesh.x,md2.mesh.y)]; - spct=[spct InterpFromMeshToMesh2d(md1.mesh.elements,md1.mesh.x,md1.mesh.y,thickness,md2.mesh.x,md2.mesh.y)]; - steps=[steps t*md1.timestepping.time_step]; - end - count=count+1; - end - - %As long as there are recorded time steps, spc the boundaries with velocities - if nsteps > 0 - md2.diagnostic.spcvx=md2.diagnostic.spcvx*ones(1,size(spcx,2)); - md2.diagnostic.spcvy=md2.diagnostic.spcvy*ones(1,size(spcy,2)); - md2.diagnostic.spcvz=md2.diagnostic.spcvz*ones(1,size(spcx,2)); - md2.prognostic.spcthickness=md2.prognostic.spcthickness*ones(1,size(spct,2)); - md2.diagnostic.spcvx(find(md2.mesh.vertexonboundary),:)=spcx(find(md2.mesh.vertexonboundary),:); - md2.diagnostic.spcvy(find(md2.mesh.vertexonboundary),:)=spcy(find(md2.mesh.vertexonboundary),:); - md2.diagnostic.spcvz(find(md2.mesh.vertexonboundary),:)=0; - md2.prognostic.spcthickness(find(md2.mesh.vertexonboundary),:)=spct(find(md2.mesh.vertexonboundary),:); - md2.diagnostic.spcvx=[md2.diagnostic.spcvx; steps]; - md2.diagnostic.spcvy=[md2.diagnostic.spcvy; steps]; - md2.diagnostic.spcvz=[md2.diagnostic.spcvz; steps]; - md2.prognostic.spcthickness=[md2.prognostic.spcthickness; steps]; - end - - %Diagnostic. Don't spc the icefront vertices. - if ~isnan(md2.diagnostic.icefront) - md1s=modelextract(md1,area); - %md2.diagnostic.icefront=[md2.mesh.segments 2]; - e2=md2.mesh.segments(:,end); - e1=md1s.mesh.segments(:,end); - - pload = nan*ones(size(md1s.mesh.elements,1),1); - pload(md1s.diagnostic.icefront(:,end-1))=md1s.diagnostic.icefront(:,end); - - x2=mean(md2.mesh.x(md2.mesh.elements(e2,:)),2); - y2=mean(md2.mesh.y(md2.mesh.elements(e2,:)),2); - x1=mean(md1s.mesh.x(md1s.mesh.elements),2); - y1=mean(md1s.mesh.y(md1s.mesh.elements),2); - - pload2=griddata(x1,y1,pload,x2,y2,'nearest'); - md2.diagnostic.icefront=[md2.mesh.segments(~isnan(pload2),:) pload2(~isnan(pload2))]; - md2.diagnostic.spcvx(unique(md2.diagnostic.icefront(:,1:2)),:)=nan; - md2.diagnostic.spcvy(unique(md2.diagnostic.icefront(:,1:2)),:)=nan; - md2.diagnostic.spcvz(unique(md2.diagnostic.icefront(:,1:2)),:)=nan; - md2.prognostic.spcthickness(unique(md2.diagnostic.icefront(:,1:2)),:)=nan; - end - - %Clear results fields - if isstruct(md1.results), - md2.results=[]; - end - Index: sm/trunk-jpl/src/m/model/shear2d.m =================================================================== --- /issm/trunk-jpl/src/m/model/shear2d.m (revision 13007) +++ (revision ) @@ -1,23 +1,0 @@ -function [sx,sy,sxy,s]=shear2d(md) -%SHEAR2D - computes 2d strain rate -% -% This routine computes the strain rate of 2d models -% -% Usage: -% [sx,sy,sxy,s]=shear2d(md); -% s=shear2d(md); - -[alpha beta]=GetNodalFunctionsCoeff(md.mesh.elements,md.mesh.x,md.mesh.y); - -summation=[1;1;1]; -sx=(md.initialization.vx(md.mesh.elements).*alpha)*summation; -uy=(md.initialization.vx(md.mesh.elements).*beta)*summation; -vx=(md.initialization.vy(md.mesh.elements).*alpha)*summation; -sy=(md.initialization.vy(md.mesh.elements).*beta)*summation; -sxy=(uy+vx)/2; -s=sqrt(sx.^2+sy.^2+sxy.^2+sx.*sy); - -%if user requested only one output, it must be the norm -if nargout==1, - sx=s; -end Index: sm/trunk-jpl/src/m/model/sia.m =================================================================== --- /issm/trunk-jpl/src/m/model/sia.m (revision 13007) +++ (revision ) @@ -1,26 +1,0 @@ -function [velx,vely,vel]=sia(md) -%BALVEL - computation of Shallow Ice velocities -% -% This routine uses the model of Hutter to compute the velocities -% of a 2d model using the surface slope -% -% Usage: -% [velx,vely,vel]=sia(md) - -if md.mesh.dimension~=2, - error('Only 2d meshes are allowed to compute velocity balances'); -end - -%Get slope -[sx,sy,s]=slope(md); - -%Average thickness and B over all elements. -summer=[1;1;1]; -hel=md.geometry.thickness(md.mesh.elements)*summer/3; -Bel=md.B(md.mesh.elements)*summer/3; - -Ael=Bel.^(-3); - -velx=-2*(md.materials.rho_ice*md.constants.g)^3*s.^2.*sx.*Ael/4.*hel.^4; -vely=-2*(md.materials.rho_ice*md.constants.g)^3*s.^2.*sy.*Ael/4.*hel.^4; -vel=sqrt(velx.^2+vely.^2); Index: sm/trunk-jpl/src/m/model/slope.m =================================================================== --- /issm/trunk-jpl/src/m/model/slope.m (revision 13007) +++ (revision ) @@ -1,32 +1,0 @@ -function [sx,sy,s]=slope(md) -%SLOPE - compute the surface slope -% -% Usage: -% [sx,sy,s]=slope(md) - -%load some variables (it is much faster if the variab;es are loaded from md once for all) -if (md.mesh.dimension==2), - numberofelements=md.mesh.numberofelements; - numberofnodes=md.mesh.numberofvertices; - index=md.mesh.elements; - x=md.mesh.x; y=md.mesh.y; -else - numberofelements=md.mesh.numberofelements2d; - numberofnodes=md.mesh.numberofvertices2d; - index=md.mesh.elements2d; - x=md.mesh.x2d; y=md.mesh.y2d; -end - -%compute nodal functions coefficients N(x,y)=alpha x + beta y + gamma -[alpha beta]=GetNodalFunctionsCoeff(index,x,y); - -summation=[1;1;1]; -sx=(md.geometry.surface(index).*alpha)*summation; -sy=(md.geometry.surface(index).*beta)*summation; -s=sqrt(sx.^2+sy.^2); - -if md.mesh.dimension==3, - sx=project3d(md,'vector',sx,'type','element'); - sy=project3d(md,'vector',sy,'type','element'); - s=sqrt(sx.^2+sy.^2); -end Index: sm/trunk-jpl/src/m/model/thicknessevolution.m =================================================================== --- /issm/trunk-jpl/src/m/model/thicknessevolution.m (revision 13007) +++ (revision ) @@ -1,28 +1,0 @@ -function dhdt=thicknessevolution(md) -%THICKNESSEVOLUTION - compute the new thickness of a model after ∆t -% -% This routine compute the new thickness of a model after a time step -% according to the following formula: -% dh/dt=-div(Hu) -% -% Usage: -% dhdt=thicknessevolution(md) - -if (length(md.initialization.vx)~=md.mesh.numberofvertices)|(length(md.initialization.vy)~=md.mesh.numberofvertices) - error(['thicknessevolution error message: vx and vy should have a length of ' num2str(md.mesh.numberofvertices)]) -end - -%load some variables -H=md.geometry.thickness; -vx=md.initialization.vx; -vy=md.initialization.vy; -index=md.mesh.elements; - -%compute nodal functions coefficients N(x,y)=alpha x + beta y + gamma -[alpha beta]=GetNodalFunctionsCoeff(md.mesh.elements,md.mesh.x,md.mesh.y); - -%compute dhdt=div(Hu) -summation=1/3*ones(3,1); -dhdt=(vx(index)*summation).*sum( H(index).*alpha,2) + (vy(index)*summation).*sum(H(index).*beta,2) ... - + ( H(index)*summation).*sum(vx(index).*alpha,2) + ( H(index)*summation).*sum(vy(index).*beta,2); -dhdt=-dhdt; Index: sm/trunk-jpl/src/m/model/tres.m =================================================================== --- /issm/trunk-jpl/src/m/model/tres.m (revision 13007) +++ (revision ) @@ -1,119 +1,0 @@ -function md=tres(md,string) -%TRES - transfer results results to corresponding model fields. -% -% Usage: md=tres(md,string) -% -% Example: md=tres(md,'diagnostic'); - -%check number of arguments - -if strcmpi(string,'diagnostic'), - if md.mesh.dimension==2, - md.initialization.vx=md.results.DiagnosticSolution.Vx; - md.initialization.vy=md.results.DiagnosticSolution.Vy; - else - md.initialization.vx=md.results.DiagnosticSolution.Vx; - md.initialization.vy=md.results.DiagnosticSolution.Vy; - md.initialization.vz=md.results.DiagnosticSolution.Vz; - end - md.initialization.vel=md.results.DiagnosticSolution.Vel; - - if isfield(md.results.DiagnosticSolution,'Pressure'), - md.initialization.pressure=md.results.DiagnosticSolution.Pressure; - end - if md.rifts.numrifts, - if isfield(md.results.DiagnosticSolution,'riftproperties'), - md.rifts.riftproperties=md.results.DiagnosticSolution.riftproperties; - end - end - if md.inversion.iscontrol==1, - for control_parameters=md.inversion.control_parameters - %Will need to be updated... good luck ;) - md.(EnumToModelField(control_parameters))=md.results.DiagnosticSolution.(EnumToString(control_parameters)); - end - end - -elseif strcmpi(string,'dakota'), - md.qmu.results=md.results.dakota; - -elseif strcmpi(string,'flaim'), - md.flaim.solution=md.results.FlaimSolution.solution; - md.flaim.quality =md.results.FlaimSolution.quality; - -elseif strcmpi(string,'transient'), - results=md.results.TransientSolution; - results2.Vel=NaN; - count=1; - for i=1:length(results), - if ~isempty(md.results.TransientSolution(i).Vel), - results2(count).Vel=md.results.TransientSolution(i).Vel; - results2(count).Surface=md.results.TransientSolution(i).Surface; - results2(count).Thickness=md.results.TransientSolution(i).Thickness; - results2(count).Bed=md.results.TransientSolution(i).Bed; - results2(count).Vx=md.results.TransientSolution(i).Vx; - results2(count).Vy=md.results.TransientSolution(i).Vy; - results2(count).time=md.results.TransientSolution(i).time; - results2(count).step=md.results.TransientSolution(i).step; - if ~strcmpi(md.groundingline.migration,'None'), - results2(count).ElementOnIceShelf=md.results.TransientSolution(i).ElementOnIceShelf; - end - count=count+1; - end - end - md.results.TransientSolution=results2; - clear results,results2; -elseif strcmpi(string,'steadystate'), - md.initialization.vx=md.results.SteadystateSolution.Vx; - md.initialization.vy=md.results.SteadystateSolution.Vy; - if isfield(md.results.SteadystateSolution,'Vz'), - md.initialization.vz=md.results.SteadystateSolution.Vz; - end - - md.initialization.vel=md.results.SteadystateSolution.Vel; - md.initialization.pressure=md.results.SteadystateSolution.Pressure; - md.initialization.temperature=md.results.SteadystateSolution.Temperature; - md.basalforcings.melting_rate=md.results.SteadystateSolution.BasalforcingsMeltingRate; - - if md.inversion.iscontrol==1, - for control_parameters=md.inversion.control_parameters - md.(EnumToModelField(control_parameters))=md.results.SteadystateSolution.(EnumToString(control_parameters)); - end - end - -elseif strcmpi(string,'thermal'), - md.initialization.temperature=md.results.ThermalSolution.Temperature; - md.basalforcings.melting_rate=md.results.ThermalSolution.BasalMeltingRate; -elseif strcmpi(string,'hydrology'), - md.initialization.watercolumn=md.results.HydrologySolution.Watercolumn; - -else - error(['tres error message: analysis ' string ' not supported yet!']); -end -end -function string=EnumToModelField(enum) % {{{ - %ENUMTOMODELFIELD - output string of model field associated to enum - % - % Usage: - % string=EnumToModelField(enum) - - disp('Warning: EnumToModelField is deprecated, it cannot work with new model definition. This function will be removed in the future'); - - switch enum, - - case ThicknessEnum(), string='thickness'; return - case FrictionCoefficientEnum(), string='drag_coefficient'; return - case MaterialsRheologyBEnum(), string='rheology_B'; return - case MaterialsRheologyBbarEnum(), string='rheology_B'; return - case BalancethicknessThickeningRateEnum: string='dhdt'; return - case VxEnum(), string='vx'; return - case InversionVxObsEnum(), string='vx_obs'; return - case VyEnum(), string='vy'; return - case InversionVyObsEnum(), string='vy_obs'; return - case BasalforcingsMeltingRateEnum(), string='basal_melting_rate'; return - case SurfaceforcingsAccumulationRateEnum(), string='surface_accumulation_rate'; return - case SurfaceforcingsAblationRateEnum(), string='surface_ablation_rate'; return - case SurfaceforcingsMassBalanceEnum(), string='surface_mass_balance'; return - otherwise, error(['Enum ' num2str(enum) ' not found associated to any model field']); - - end - end % }}}