Index: /issm/trunk-jpl/src/m/contrib/massbalance/contourmassbalance.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/massbalance/contourmassbalance.m (revision 13006) +++ /issm/trunk-jpl/src/m/contrib/massbalance/contourmassbalance.m (revision 13006) @@ -0,0 +1,45 @@ +function dhdt=contourmassbalance(md,file) +%CONTOURMASSBALANCE - compute the mass balance on a contour +% +% Usage: +% dhdt=contourmassbalance(md,file) + +%some checks +if nargin~=2, + help contourmassbalance + error('contourmassbalance error message: bad usage'); +end +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 +if ~exist(file), + error(['thicknessevolution error message: file ' file ' not found']); +end + +%Get segments enveloping contour +segments=contourenvelope(md,file); +%md.diagnostic.icefront=segments; plotmodel(md,'data','pressureload','expdisp',file); + +%get flag list of elements and nodes inside the contour +nodein=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,file,'node',1); +elemin=(sum(nodein(md.mesh.elements),2)==size(md.mesh.elements,2)); + +%conputing Mass flux +x=md.mesh.x; +y=md.mesh.y; +vx=mean(md.initialization.vx(segments(:,1:end-1)),2); +vy=mean(md.initialization.vy(segments(:,1:end-1)),2); +H=mean(md.geometry.thickness(segments(:,1:end-1)),2); +nx=cos(atan2((x(segments(:,1))-x(segments(:,2))) , (y(segments(:,2))-y(segments(:,1))))); +ny=sin(atan2((x(segments(:,1))-x(segments(:,2))) , (y(segments(:,2))-y(segments(:,1))))); +L=sqrt((x(segments(:,1))-x(segments(:,2))).^2+(y(segments(:,2))-y(segments(:,1))).^2); +flux = - md.materials.rho_ice*sum(L.*H.*(vx.*nx+vy.*ny)); %outflux is negative! +disp(['mass outflux on ' file ' = ' num2str(-flux/10^9) ' Gt/yr']); +areas=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y); +dhdt=flux/(sum(areas(find(elemin)))*md.materials.rho_ice); +disp(['dhdt on ' file ' (Flux method) = ' num2str(dhdt) ' m/yr']); + +dhdt=thicknessevolution(md); +in=find(elemin); +dhdt=sum(dhdt(in).*areas(in))/sum(areas(in)); +disp(['dhdt on ' file ' (divHV method) = ' num2str(dhdt) ' m/yr']); Index: sm/trunk-jpl/src/m/model/contourenvelope.m =================================================================== --- /issm/trunk-jpl/src/m/model/contourenvelope.m (revision 13005) +++ (revision ) @@ -1,139 +1,0 @@ -function segments=contourenvelope(md,varargin) -%CONTOURENVELOPE - build a set of segments enveloping a contour .exp -% -% Usage: -% segments=contourenvelope(md,varargin) -% -% Example: -% segments=contourenvelope(md,'Stream.exp'); -% segments=contourenvelope(md,md.mask.elementonfloatingice) -% segments=contourenvelope(md); - -%some checks -if nargin>2, - help contourenvelope - error('contourenvelope error message: bad usage'); -end -if nargin==2, - flags=varargin{1}; - - if ischar(flags), - file=flags; - file=varargin{1}; - if ~exist(file), - error(['contourenvelope error message: file ' file ' not found']); - end - isfile=1; - elseif isnumeric(flags), - %do nothing for now - isfile=0; - else - error('contourenvelope error message: second argument should a file or an elements flag'); - end -end - -%Now, build the connectivity tables for this mesh. -%Computing connectivity -if (size(md.mesh.vertexconnectivity,1)~=md.mesh.numberofvertices & size(md.mesh.vertexconnectivity,1)~=md.mesh.numberofvertices2d), - md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices); -end -if (size(md.mesh.elementconnectivity,1)~=md.mesh.numberofelements & size(md.mesh.elementconnectivity,1)~=md.mesh.numberofelements2d), - md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity); -end - -%get nodes inside profile -mesh.elementconnectivity=md.mesh.elementconnectivity; -if md.mesh.dimension==2; - mesh.elements=md.mesh.elements; - mesh.x=md.mesh.x; - mesh.y=md.mesh.y; - mesh.numberofvertices=md.mesh.numberofvertices; - mesh.numberofelements=md.mesh.numberofelements; -else - mesh.elements=md.mesh.elements2d; - mesh.x=md.mesh.x2d; - mesh.y=md.mesh.y2d; - mesh.numberofvertices=md.mesh.numberofvertices2d; - mesh.numberofelements=md.mesh.numberofelements2d; -end - -if nargin==2, - - if isfile, - %get flag list of elements and nodes inside the contour - nodein=ContourToMesh(mesh.elements,mesh.x,mesh.y,file,'node',1); - elemin=(sum(nodein(mesh.elements),2)==size(mesh.elements,2)); - %modify element connectivity - elemout=find(~elemin); - mesh.elementconnectivity(elemout,:)=0; - mesh.elementconnectivity(find(ismember(mesh.elementconnectivity,elemout)))=0; - else - %get flag list of elements and nodes inside the contour - nodein=zeros(mesh.numberofvertices,1); - elemin=zeros(mesh.numberofelements,1); - - pos=find(flags); - elemin(pos)=1; - nodein(mesh.elements(pos,:))=1; - - %modify element connectivity - elemout=find(~elemin); - mesh.elementconnectivity(elemout,:)=0; - mesh.elementconnectivity(find(ismember(mesh.elementconnectivity,elemout)))=0; - end -end - -%Find element on boundary -%First: find elements on the boundary of the domain -flag=mesh.elementconnectivity; -if nargin==2, - flag(find(flag))=elemin(flag(find(flag))); -end -elementonboundary=double(prod(flag,2)==0 & sum(flag,2)>0); - -%Find segments on boundary -pos=find(elementonboundary); -num_segments=length(pos); -segments=zeros(num_segments,3); -count=1; - -for i=1:num_segments, - el1=pos(i); - els2=mesh.elementconnectivity(el1,find(mesh.elementconnectivity(el1,:))); - if length(els2)>1, - flag=intersect(mesh.elements(els2(1),:),mesh.elements(els2(2),:)); - nods1=mesh.elements(el1,:); - nods1(find(nods1==flag))=[]; - segments(count,:)=[nods1 el1]; - - ord1=find(nods1(1)==mesh.elements(el1,:)); - ord2=find(nods1(2)==mesh.elements(el1,:)); - - %swap segment nodes if necessary - if ( (ord1==1 & ord2==2) | (ord1==2 & ord2==3) | (ord1==3 & ord2==1) ), - temp=segments(count,1); - segments(count,1)=segments(count,2); - segments(count,2)=temp; - end - segments(count,1:2)=fliplr(segments(count,1:2)); - count=count+1; - else - nods1=mesh.elements(el1,:); - flag=setdiff(nods1,mesh.elements(els2,:)); - for j=1:3, - nods=nods1; nods(j)=[]; - if any(ismember(flag,nods)), - segments(count,:)=[nods el1]; - ord1=find(nods(1)==mesh.elements(el1,:)); - ord2=find(nods(2)==mesh.elements(el1,:)); - if ( (ord1==1 & ord2==2) | (ord1==2 & ord2==3) | (ord1==3 & ord2==1) ), - temp=segments(count,1); - segments(count,1)=segments(count,2); - segments(count,2)=temp; - end - segments(count,1:2)=fliplr(segments(count,1:2)); - count=count+1; - end - end - end -end Index: sm/trunk-jpl/src/m/model/contourmassbalance.m =================================================================== --- /issm/trunk-jpl/src/m/model/contourmassbalance.m (revision 13005) +++ (revision ) @@ -1,45 +1,0 @@ -function dhdt=contourmassbalance(md,file) -%CONTOURMASSBALANCE - compute the mass balance on a given profile -% -% Usage: -% dhdt=contourmassbalance(md,file) - -%some checks -if nargin~=2, - help contourmassbalance - error('contourmassbalance error message: bad usage'); -end -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 -if ~exist(file), - error(['thicknessevolution error message: file ' file ' not found']); -end - -%Get segments enveloping contour -segments=contourenvelope(md,file); -%md.diagnostic.icefront=segments; plotmodel(md,'data','pressureload','expdisp',file); - -%get flag list of elements and nodes inside the contour -nodein=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,file,'node',1); -elemin=(sum(nodein(md.mesh.elements),2)==size(md.mesh.elements,2)); - -%conputing Mass flux -x=md.mesh.x; -y=md.mesh.y; -vx=mean(md.initialization.vx(segments(:,1:end-1)),2); -vy=mean(md.initialization.vy(segments(:,1:end-1)),2); -H=mean(md.geometry.thickness(segments(:,1:end-1)),2); -nx=cos(atan2((x(segments(:,1))-x(segments(:,2))) , (y(segments(:,2))-y(segments(:,1))))); -ny=sin(atan2((x(segments(:,1))-x(segments(:,2))) , (y(segments(:,2))-y(segments(:,1))))); -L=sqrt((x(segments(:,1))-x(segments(:,2))).^2+(y(segments(:,2))-y(segments(:,1))).^2); -flux = - md.materials.rho_ice*sum(L.*H.*(vx.*nx+vy.*ny)); %outflux is negative! -disp(['mass outflux on ' file ' = ' num2str(-flux/10^9) ' Gt/yr']); -areas=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y); -dhdt=flux/(sum(areas(find(elemin)))*md.materials.rho_ice); -disp(['dhdt on ' file ' (Flux method) = ' num2str(dhdt) ' m/yr']); - -dhdt=thicknessevolution(md); -in=find(elemin); -dhdt=sum(dhdt(in).*areas(in))/sum(areas(in)); -disp(['dhdt on ' file ' (divHV method) = ' num2str(dhdt) ' m/yr']); Index: /issm/trunk-jpl/src/m/model/inversions/parametercontrolB.m =================================================================== --- /issm/trunk-jpl/src/m/model/inversions/parametercontrolB.m (revision 13006) +++ /issm/trunk-jpl/src/m/model/inversions/parametercontrolB.m (revision 13006) @@ -0,0 +1,122 @@ +function md=parametercontrolB(md,varargin), +%PARAMETERCONTROLB - parameterization for control method on B +% +% It is possible to specify the number of steps, values for the +% minimum and maximum values of B, the +% kind of cm_responses to use or the the optscal. +% +% Usage: +% md=parametercontrolB(md,varargin) +% +% Example: +% md=parametercontrolB(md) +% md=parametercontrolB(md,'nsteps',20,'cm_responses',0) +% md=parametercontrolB(md,'cm_min',10,'cm_max',10^8,'cm_jump',0.99,'maxiter',20) +% md=parametercontrolB(md,eps_cm',10^-4,'optscal',[10^7 10^8]) +% +% See also PARAMETERCONTROLDRAG + +%process options +options=pairoptions(varargin{:}); + +%control type +md.inversion.control_parameters={'MaterialsRheologyBbar'}; + +%weights +weights=getfieldvalue(options,'weights',ones(md.mesh.numberofvertices,1)); +if (length(weights)~=md.mesh.numberofvertices) + md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,1); +else + md.inversion.cost_functions_coefficients=weights; +end + +%nsteps +nsteps=getfieldvalue(options,'nsteps',100); +if (length(nsteps)~=1 | nsteps<=0 | floor(nsteps)~=nsteps) + md.inversion.nsteps=100; +else + md.inversion.nsteps=nsteps; +end + + +%cm_min +cm_min=getfieldvalue(options,'cm_min',paterson(273.15+5)*ones(md.mesh.numberofvertices,1)); +if (length(cm_min)==1) + md.inversion.min_parameters=cm_min*ones(md.mesh.numberofvertices,1); +elseif (length(cm_min)==md.mesh.numberofvertices) + md.inversion.min_parameters=cm_min; +else + md.inversion.min_parameters=cm_min; +end + +%cm_max +cm_max=getfieldvalue(options,'cm_max',paterson(273.15-70)*ones(md.mesh.numberofvertices,1)); +if (length(cm_max)==1) + md.inversion.max_parameters=cm_max*ones(md.mesh.numberofvertices,1); +elseif (length(cm_max)==md.mesh.numberofvertices) + md.inversion.max_parameters=cm_max; +else + md.inversion.max_parameters=cm_max; +end + +%eps_cm +eps_cm=getfieldvalue(options,'eps_cm',NaN); +if (length(eps_cm)~=1 | eps_cm<0 ) + md.inversion.cost_function_threshold=NaN; +else + md.inversion.cost_function_threshold=eps_cm; +end + +%maxiter +maxiter=getfieldvalue(options,'maxiter',10*ones(md.inversion.nsteps,1)); +if (any(maxiter<0) | any(floor(maxiter)~=maxiter)) + md.inversion.maxiter_per_step=10*ones(md.inversion.nsteps,1); +else + md.inversion.maxiter_per_step=repmat(maxiter(:),md.inversion.nsteps,1); + md.inversion.maxiter_per_step(md.inversion.nsteps+1:end)=[]; +end + +%cm_jump +cm_jump=getfieldvalue(options,'cm_jump',0.9*ones(md.inversion.nsteps,1)); +if ~isreal(cm_jump) + md.inversion.step_threshold=0.9*ones(md.inversion.nsteps,1); +else + md.inversion.step_threshold=repmat(cm_jump(:),md.inversion.nsteps,1); + md.inversion.step_threshold(md.inversion.nsteps+1:end)=[]; +end + +%cm_responses +found=0; +if exist(options,'cm_responses'), + cm_responses=getfieldvalue(options,'cm_responses'); + if ~any(~ismember(cm_responses,[ 101:105])), + md.inversion.cost_functions=repmat(cm_responses(:),md.inversion.nsteps,1); + md.inversion.cost_functions(md.inversion.nsteps+1:end)=[]; + found=1; + end +end +if ~found + third=ceil(md.inversion.nsteps/3); + md.inversion.cost_functions=[... + 103*ones(third,1);... + 101*ones(third,1);... + repmat([101;101;103;101],third,1)... + ]; + md.inversion.cost_functions(md.inversion.nsteps+1:end)=[]; +end + +%optscal +found=0; +if exist(options,'optscal'), + optscal=getfieldvalue(options,'optscal'); + if ~any(optscal<0), + md.inversion.gradient_scaling=repmat(optscal(:),md.inversion.nsteps,1); + md.inversion.gradient_scaling(md.inversion.nsteps+1:end)=[]; + found=1; + end +end +if ~found + third=ceil(md.inversion.nsteps/3); + md.inversion.gradient_scaling=[2*10^8*ones(3,1);10^8*ones(third-3,1);10^7*ones(2*third,1);]; + md.inversion.gradient_scaling(md.inversion.nsteps+1:end)=[]; +end Index: /issm/trunk-jpl/src/m/model/inversions/parametercontroldrag.m =================================================================== --- /issm/trunk-jpl/src/m/model/inversions/parametercontroldrag.m (revision 13006) +++ /issm/trunk-jpl/src/m/model/inversions/parametercontroldrag.m (revision 13006) @@ -0,0 +1,121 @@ +function md=parametercontroldrag(md,varargin), +%PARAMETERCONTROLDRAG - parameterization for control method on drag +% +% It is possible to specify the number of steps, values for the +% minimum and maximum values of the drag, the +% kind of cm_responses to use or the the optscal. +% +% Usage: +% md=parametercontroldrag(md,varargin) +% +% Example: +% md=parametercontroldrag(md) +% md=parametercontroldrag(md,'nsteps',20,'cm_responses',0) +% md=parametercontroldrag(md,'cm_min',1,'cm_max',150,'cm_jump',0.99,'maxiter',20) +% md=parametercontroldrag(md,eps_cm',10^-4,'optscal',[10^7 10^8]) +% +% See also PARAMETERCONTROLB + +%process options +options=pairoptions(varargin{:}); + +%control type +md.inversion.control_parameters={'FrictionCoefficient'}; + +%weights +weights=getfieldvalue(options,'weights',ones(md.mesh.numberofvertices,1)); +if (length(weights)~=md.mesh.numberofvertices) + md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,1); +else + md.inversion.cost_functions_coefficients=weights; +end + +%nsteps +nsteps=getfieldvalue(options,'nsteps',100); +if (length(nsteps)~=1 | nsteps<=0 | floor(nsteps)~=nsteps) + md.inversion.nsteps=100; +else + md.inversion.nsteps=nsteps; +end + +%cm_min +cm_min=getfieldvalue(options,'cm_min',1*ones(md.mesh.numberofvertices,1)); +if (length(cm_min)==1) + md.inversion.min_parameters=cm_min*ones(md.mesh.numberofvertices,1); +elseif (length(cm_min)==md.mesh.numberofvertices) + md.inversion.min_parameters=cm_min; +else + md.inversion.min_parameters=cm_min; +end + +%cm_max +cm_max=getfieldvalue(options,'cm_max',250*ones(md.mesh.numberofvertices,1)); +if (length(cm_max)==1) + md.inversion.max_parameters=cm_max*ones(md.mesh.numberofvertices,1); +elseif (length(cm_max)==md.mesh.numberofvertices) + md.inversion.max_parameters=cm_max; +else + md.inversion.max_parameters=cm_max; +end + +%eps_cm +eps_cm=getfieldvalue(options,'eps_cm',NaN); +if (length(eps_cm)~=1 | eps_cm<0 ) + md.inversion.cost_function_threshold=NaN; +else + md.inversion.cost_function_threshold=eps_cm; +end + +%maxiter +maxiter=getfieldvalue(options,'maxiter',10*ones(md.inversion.nsteps,1)); +if (any(maxiter<0) | any(floor(maxiter)~=maxiter)) + md.inversion.maxiter_per_step=10*ones(md.inversion.nsteps,1); +else + md.inversion.maxiter_per_step=repmat(maxiter(:),md.inversion.nsteps,1); + md.inversion.maxiter_per_step(md.inversion.nsteps+1:end)=[]; +end + +%cm_jump +cm_jump=getfieldvalue(options,'cm_jump',0.8*ones(md.inversion.nsteps,1)); +if ~isreal(cm_jump) + md.inversion.step_threshold=0.8*ones(md.inversion.nsteps,1); +else + md.inversion.step_threshold=repmat(cm_jump(:),md.inversion.nsteps,1); + md.inversion.step_threshold(md.inversion.nsteps+1:end)=[]; +end + +%cm_responses +found=0; +if exist(options,'cm_responses'), + cm_responses=getfieldvalue(options,'cm_responses'); + if ~any(~ismember(cm_responses,[101 105])) + md.inversion.cost_functions=repmat(cm_responses(:),md.inversion.nsteps,1); + md.inversion.cost_functions(md.inversion.nsteps+1:end)=[]; + found=1; + end +end +if ~found + third=ceil(md.inversion.nsteps/3); + md.inversion.cost_functions=[... + 103*ones(third,1);... + 101*ones(third,1);... + repmat([101;101;103;101],third,1)... + ]; + md.inversion.cost_functions(md.inversion.nsteps+1:end)=[]; +end + +%optscal +found=0; +if exist(options,'optscal'), + optscal=getfieldvalue(options,'optscal'); + if ~any(optscal<0), + md.inversion.gradient_scaling=repmat(optscal(:),md.inversion.nsteps,1); + md.inversion.gradient_scaling(md.inversion.nsteps+1:end)=[]; + found=1; + end +end +if ~found + third=ceil(md.inversion.nsteps/3); + md.inversion.gradient_scaling=[50*ones(3,1);15*ones(third-3,1);10*ones(third,1);repmat([10;10;20;10],third,1)]; + md.inversion.gradient_scaling(md.inversion.nsteps+1:end)=[]; +end Index: sm/trunk-jpl/src/m/model/parametercontroloptimization.m =================================================================== --- /issm/trunk-jpl/src/m/model/parametercontroloptimization.m (revision 13005) +++ (revision ) @@ -1,70 +1,0 @@ -function md=parametercontroloptimization(md,varargin), -%PARAMETERCONTROLOPTIMIZATION - parameterization for control method on drag -% -% It is possible to specify the number of steps, values for the -% minimum and maximum values of the drag, specify the -% kind of fit to use or the the optscal. -% -% Usage: -% md=parametercontroloptimization(md,varargin) -% -% Example: -% md=parametercontroloptimization(md,'nsteps',6) - -%process options -options=pairoptions(varargin{:}); - -%Copy model -md2=md; - -%Hard coded parameters -%variable -cmax_max=1000; -cm_maxs=linspace(cmax_max/10,cmax_max,3); -optscal_max=1000; -optscals=linspace(optscal_max/10,optscal_max,3); -fits=[0 2]; -%Kept constant -md2.verbose=0; -md2.cm_min=1; -md2.eps_cm=NaN; -md2.nsteps=getfieldvalue(options,'nsteps',5); -md2.inversion.control_parameters=getfieldvalue(options,'md2.inversion.control_parameters',{'FrictionCoefficient'}); -md2.maxiter=10*ones(md2.nsteps,1); -md2.cm_jump=0.99*ones(md2.nsteps,1); -md2.inversion.iscontrol=1; -md2.weights=ones(md2.numberofnodes,1); - -%loop over the set of parameters -best=0; -count=1; -for fit=fits, - md2.fit=fit*ones(md2.nsteps,1); - for cm_max=cm_maxs; - md2.cm_max=cm_max; - for optscal=optscals; - md2.optscal=optscal*ones(md2.nsteps,1); - - %current run - disp(sprintf('\n Step %i/%i, fit=%i, cm_max=%g, optscal=%g\n',count,length(cm_maxs)*length(optscals)*length(fits),fit,cm_max,optscal)) - md2=solve(md2,'analysis_type','DiagnosticAnalysis'); - count=count+1; - - %Check misfit difference - rel=(md2.results.DiagnosticAnalysis.J(1)-md2.results.DiagnosticAnalysis.J(end))/md2.results.DiagnosticAnalysis.J(1); - disp([' ΔJ/J=' num2str(rel*100) '% for fit=' num2str(fit) ',cm_max=' num2str(cm_max) ', and optscal=' num2str(optscal) ]); - if (rel>best), - disp(sprintf(' --> current best')) - fit_final =md2.fit; - cmmax_final =md2.cm_max; - optscal_final=md2.optscal; - best=rel; - end - end - end -end - -%load final parameters onto initial model -md.fit =fit_final; -md.inversion.max_parameters =cmmax_final; -md.inversion.gradient_scaling=optscal_final; Index: /issm/trunk-jpl/src/m/model/parameterization/contourenvelope.m =================================================================== --- /issm/trunk-jpl/src/m/model/parameterization/contourenvelope.m (revision 13006) +++ /issm/trunk-jpl/src/m/model/parameterization/contourenvelope.m (revision 13006) @@ -0,0 +1,139 @@ +function segments=contourenvelope(md,varargin) +%CONTOURENVELOPE - build a set of segments enveloping a contour .exp +% +% Usage: +% segments=contourenvelope(md,varargin) +% +% Example: +% segments=contourenvelope(md,'Stream.exp'); +% segments=contourenvelope(md,md.mask.elementonfloatingice) +% segments=contourenvelope(md); + +%some checks +if nargin>2, + help contourenvelope + error('contourenvelope error message: bad usage'); +end +if nargin==2, + flags=varargin{1}; + + if ischar(flags), + file=flags; + file=varargin{1}; + if ~exist(file), + error(['contourenvelope error message: file ' file ' not found']); + end + isfile=1; + elseif isnumeric(flags), + %do nothing for now + isfile=0; + else + error('contourenvelope error message: second argument should a file or an elements flag'); + end +end + +%Now, build the connectivity tables for this mesh. +%Computing connectivity +if (size(md.mesh.vertexconnectivity,1)~=md.mesh.numberofvertices & size(md.mesh.vertexconnectivity,1)~=md.mesh.numberofvertices2d), + md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices); +end +if (size(md.mesh.elementconnectivity,1)~=md.mesh.numberofelements & size(md.mesh.elementconnectivity,1)~=md.mesh.numberofelements2d), + md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity); +end + +%get nodes inside profile +mesh.elementconnectivity=md.mesh.elementconnectivity; +if md.mesh.dimension==2; + mesh.elements=md.mesh.elements; + mesh.x=md.mesh.x; + mesh.y=md.mesh.y; + mesh.numberofvertices=md.mesh.numberofvertices; + mesh.numberofelements=md.mesh.numberofelements; +else + mesh.elements=md.mesh.elements2d; + mesh.x=md.mesh.x2d; + mesh.y=md.mesh.y2d; + mesh.numberofvertices=md.mesh.numberofvertices2d; + mesh.numberofelements=md.mesh.numberofelements2d; +end + +if nargin==2, + + if isfile, + %get flag list of elements and nodes inside the contour + nodein=ContourToMesh(mesh.elements,mesh.x,mesh.y,file,'node',1); + elemin=(sum(nodein(mesh.elements),2)==size(mesh.elements,2)); + %modify element connectivity + elemout=find(~elemin); + mesh.elementconnectivity(elemout,:)=0; + mesh.elementconnectivity(find(ismember(mesh.elementconnectivity,elemout)))=0; + else + %get flag list of elements and nodes inside the contour + nodein=zeros(mesh.numberofvertices,1); + elemin=zeros(mesh.numberofelements,1); + + pos=find(flags); + elemin(pos)=1; + nodein(mesh.elements(pos,:))=1; + + %modify element connectivity + elemout=find(~elemin); + mesh.elementconnectivity(elemout,:)=0; + mesh.elementconnectivity(find(ismember(mesh.elementconnectivity,elemout)))=0; + end +end + +%Find element on boundary +%First: find elements on the boundary of the domain +flag=mesh.elementconnectivity; +if nargin==2, + flag(find(flag))=elemin(flag(find(flag))); +end +elementonboundary=double(prod(flag,2)==0 & sum(flag,2)>0); + +%Find segments on boundary +pos=find(elementonboundary); +num_segments=length(pos); +segments=zeros(num_segments,3); +count=1; + +for i=1:num_segments, + el1=pos(i); + els2=mesh.elementconnectivity(el1,find(mesh.elementconnectivity(el1,:))); + if length(els2)>1, + flag=intersect(mesh.elements(els2(1),:),mesh.elements(els2(2),:)); + nods1=mesh.elements(el1,:); + nods1(find(nods1==flag))=[]; + segments(count,:)=[nods1 el1]; + + ord1=find(nods1(1)==mesh.elements(el1,:)); + ord2=find(nods1(2)==mesh.elements(el1,:)); + + %swap segment nodes if necessary + if ( (ord1==1 & ord2==2) | (ord1==2 & ord2==3) | (ord1==3 & ord2==1) ), + temp=segments(count,1); + segments(count,1)=segments(count,2); + segments(count,2)=temp; + end + segments(count,1:2)=fliplr(segments(count,1:2)); + count=count+1; + else + nods1=mesh.elements(el1,:); + flag=setdiff(nods1,mesh.elements(els2,:)); + for j=1:3, + nods=nods1; nods(j)=[]; + if any(ismember(flag,nods)), + segments(count,:)=[nods el1]; + ord1=find(nods(1)==mesh.elements(el1,:)); + ord2=find(nods(2)==mesh.elements(el1,:)); + if ( (ord1==1 & ord2==2) | (ord1==2 & ord2==3) | (ord1==3 & ord2==1) ), + temp=segments(count,1); + segments(count,1)=segments(count,2); + segments(count,2)=temp; + end + segments(count,1:2)=fliplr(segments(count,1:2)); + count=count+1; + end + end + end +end Index: sm/trunk-jpl/src/m/model/parameterization/parametercontrolB.m =================================================================== --- /issm/trunk-jpl/src/m/model/parameterization/parametercontrolB.m (revision 13005) +++ (revision ) @@ -1,122 +1,0 @@ -function md=parametercontrolB(md,varargin), -%PARAMETERCONTROLB - parameterization for control method on B -% -% It is possible to specify the number of steps, values for the -% minimum and maximum values of B, the -% kind of cm_responses to use or the the optscal. -% -% Usage: -% md=parametercontrolB(md,varargin) -% -% Example: -% md=parametercontrolB(md) -% md=parametercontrolB(md,'nsteps',20,'cm_responses',0) -% md=parametercontrolB(md,'cm_min',10,'cm_max',10^8,'cm_jump',0.99,'maxiter',20) -% md=parametercontrolB(md,eps_cm',10^-4,'optscal',[10^7 10^8]) -% -% See also PARAMETERCONTROLDRAG - -%process options -options=pairoptions(varargin{:}); - -%control type -md.inversion.control_parameters={'MaterialsRheologyBbar'}; - -%weights -weights=getfieldvalue(options,'weights',ones(md.mesh.numberofvertices,1)); -if (length(weights)~=md.mesh.numberofvertices) - md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,1); -else - md.inversion.cost_functions_coefficients=weights; -end - -%nsteps -nsteps=getfieldvalue(options,'nsteps',100); -if (length(nsteps)~=1 | nsteps<=0 | floor(nsteps)~=nsteps) - md.inversion.nsteps=100; -else - md.inversion.nsteps=nsteps; -end - - -%cm_min -cm_min=getfieldvalue(options,'cm_min',paterson(273.15+5)*ones(md.mesh.numberofvertices,1)); -if (length(cm_min)==1) - md.inversion.min_parameters=cm_min*ones(md.mesh.numberofvertices,1); -elseif (length(cm_min)==md.mesh.numberofvertices) - md.inversion.min_parameters=cm_min; -else - md.inversion.min_parameters=cm_min; -end - -%cm_max -cm_max=getfieldvalue(options,'cm_max',paterson(273.15-70)*ones(md.mesh.numberofvertices,1)); -if (length(cm_max)==1) - md.inversion.max_parameters=cm_max*ones(md.mesh.numberofvertices,1); -elseif (length(cm_max)==md.mesh.numberofvertices) - md.inversion.max_parameters=cm_max; -else - md.inversion.max_parameters=cm_max; -end - -%eps_cm -eps_cm=getfieldvalue(options,'eps_cm',NaN); -if (length(eps_cm)~=1 | eps_cm<0 ) - md.inversion.cost_function_threshold=NaN; -else - md.inversion.cost_function_threshold=eps_cm; -end - -%maxiter -maxiter=getfieldvalue(options,'maxiter',10*ones(md.inversion.nsteps,1)); -if (any(maxiter<0) | any(floor(maxiter)~=maxiter)) - md.inversion.maxiter_per_step=10*ones(md.inversion.nsteps,1); -else - md.inversion.maxiter_per_step=repmat(maxiter(:),md.inversion.nsteps,1); - md.inversion.maxiter_per_step(md.inversion.nsteps+1:end)=[]; -end - -%cm_jump -cm_jump=getfieldvalue(options,'cm_jump',0.9*ones(md.inversion.nsteps,1)); -if ~isreal(cm_jump) - md.inversion.step_threshold=0.9*ones(md.inversion.nsteps,1); -else - md.inversion.step_threshold=repmat(cm_jump(:),md.inversion.nsteps,1); - md.inversion.step_threshold(md.inversion.nsteps+1:end)=[]; -end - -%cm_responses -found=0; -if exist(options,'cm_responses'), - cm_responses=getfieldvalue(options,'cm_responses'); - if ~any(~ismember(cm_responses,[ 101:105])), - md.inversion.cost_functions=repmat(cm_responses(:),md.inversion.nsteps,1); - md.inversion.cost_functions(md.inversion.nsteps+1:end)=[]; - found=1; - end -end -if ~found - third=ceil(md.inversion.nsteps/3); - md.inversion.cost_functions=[... - 103*ones(third,1);... - 101*ones(third,1);... - repmat([101;101;103;101],third,1)... - ]; - md.inversion.cost_functions(md.inversion.nsteps+1:end)=[]; -end - -%optscal -found=0; -if exist(options,'optscal'), - optscal=getfieldvalue(options,'optscal'); - if ~any(optscal<0), - md.inversion.gradient_scaling=repmat(optscal(:),md.inversion.nsteps,1); - md.inversion.gradient_scaling(md.inversion.nsteps+1:end)=[]; - found=1; - end -end -if ~found - third=ceil(md.inversion.nsteps/3); - md.inversion.gradient_scaling=[2*10^8*ones(3,1);10^8*ones(third-3,1);10^7*ones(2*third,1);]; - md.inversion.gradient_scaling(md.inversion.nsteps+1:end)=[]; -end Index: sm/trunk-jpl/src/m/model/parameterization/parametercontroldrag.m =================================================================== --- /issm/trunk-jpl/src/m/model/parameterization/parametercontroldrag.m (revision 13005) +++ (revision ) @@ -1,121 +1,0 @@ -function md=parametercontroldrag(md,varargin), -%PARAMETERCONTROLDRAG - parameterization for control method on drag -% -% It is possible to specify the number of steps, values for the -% minimum and maximum values of the drag, the -% kind of cm_responses to use or the the optscal. -% -% Usage: -% md=parametercontroldrag(md,varargin) -% -% Example: -% md=parametercontroldrag(md) -% md=parametercontroldrag(md,'nsteps',20,'cm_responses',0) -% md=parametercontroldrag(md,'cm_min',1,'cm_max',150,'cm_jump',0.99,'maxiter',20) -% md=parametercontroldrag(md,eps_cm',10^-4,'optscal',[10^7 10^8]) -% -% See also PARAMETERCONTROLB - -%process options -options=pairoptions(varargin{:}); - -%control type -md.inversion.control_parameters={'FrictionCoefficient'}; - -%weights -weights=getfieldvalue(options,'weights',ones(md.mesh.numberofvertices,1)); -if (length(weights)~=md.mesh.numberofvertices) - md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,1); -else - md.inversion.cost_functions_coefficients=weights; -end - -%nsteps -nsteps=getfieldvalue(options,'nsteps',100); -if (length(nsteps)~=1 | nsteps<=0 | floor(nsteps)~=nsteps) - md.inversion.nsteps=100; -else - md.inversion.nsteps=nsteps; -end - -%cm_min -cm_min=getfieldvalue(options,'cm_min',1*ones(md.mesh.numberofvertices,1)); -if (length(cm_min)==1) - md.inversion.min_parameters=cm_min*ones(md.mesh.numberofvertices,1); -elseif (length(cm_min)==md.mesh.numberofvertices) - md.inversion.min_parameters=cm_min; -else - md.inversion.min_parameters=cm_min; -end - -%cm_max -cm_max=getfieldvalue(options,'cm_max',250*ones(md.mesh.numberofvertices,1)); -if (length(cm_max)==1) - md.inversion.max_parameters=cm_max*ones(md.mesh.numberofvertices,1); -elseif (length(cm_max)==md.mesh.numberofvertices) - md.inversion.max_parameters=cm_max; -else - md.inversion.max_parameters=cm_max; -end - -%eps_cm -eps_cm=getfieldvalue(options,'eps_cm',NaN); -if (length(eps_cm)~=1 | eps_cm<0 ) - md.inversion.cost_function_threshold=NaN; -else - md.inversion.cost_function_threshold=eps_cm; -end - -%maxiter -maxiter=getfieldvalue(options,'maxiter',10*ones(md.inversion.nsteps,1)); -if (any(maxiter<0) | any(floor(maxiter)~=maxiter)) - md.inversion.maxiter_per_step=10*ones(md.inversion.nsteps,1); -else - md.inversion.maxiter_per_step=repmat(maxiter(:),md.inversion.nsteps,1); - md.inversion.maxiter_per_step(md.inversion.nsteps+1:end)=[]; -end - -%cm_jump -cm_jump=getfieldvalue(options,'cm_jump',0.8*ones(md.inversion.nsteps,1)); -if ~isreal(cm_jump) - md.inversion.step_threshold=0.8*ones(md.inversion.nsteps,1); -else - md.inversion.step_threshold=repmat(cm_jump(:),md.inversion.nsteps,1); - md.inversion.step_threshold(md.inversion.nsteps+1:end)=[]; -end - -%cm_responses -found=0; -if exist(options,'cm_responses'), - cm_responses=getfieldvalue(options,'cm_responses'); - if ~any(~ismember(cm_responses,[101 105])) - md.inversion.cost_functions=repmat(cm_responses(:),md.inversion.nsteps,1); - md.inversion.cost_functions(md.inversion.nsteps+1:end)=[]; - found=1; - end -end -if ~found - third=ceil(md.inversion.nsteps/3); - md.inversion.cost_functions=[... - 103*ones(third,1);... - 101*ones(third,1);... - repmat([101;101;103;101],third,1)... - ]; - md.inversion.cost_functions(md.inversion.nsteps+1:end)=[]; -end - -%optscal -found=0; -if exist(options,'optscal'), - optscal=getfieldvalue(options,'optscal'); - if ~any(optscal<0), - md.inversion.gradient_scaling=repmat(optscal(:),md.inversion.nsteps,1); - md.inversion.gradient_scaling(md.inversion.nsteps+1:end)=[]; - found=1; - end -end -if ~found - third=ceil(md.inversion.nsteps/3); - md.inversion.gradient_scaling=[50*ones(3,1);15*ones(third-3,1);10*ones(third,1);repmat([10;10;20;10],third,1)]; - md.inversion.gradient_scaling(md.inversion.nsteps+1:end)=[]; -end Index: /issm/trunk-jpl/src/m/model/parameterization/setflowequation.m =================================================================== --- /issm/trunk-jpl/src/m/model/parameterization/setflowequation.m (revision 13006) +++ /issm/trunk-jpl/src/m/model/parameterization/setflowequation.m (revision 13006) @@ -0,0 +1,281 @@ +function md=setflowequation(md,varargin) +%SETELEMENTSTYPE - associate a solution type to each element +% +% This routine works like plotmodel: it works with an even number of inputs +% 'hutter','macayeal','pattyn','stokes' and 'fill' are the possible options +% that must be followed by the corresponding exp file or flags list +% It can either be a domain file (argus type, .exp extension), or an array of element flags. +% If user wants every element outside the domain to be +% setflowequationd, add '~' to the name of the domain file (ex: '~Pattyn.exp'); +% an empty string '' will be considered as an empty domain +% a string 'all' will be considered as the entire domain +% You can specify the type of coupling, 'penalties' or 'tiling', to use with the input 'coupling' +% +% Usage: +% md=setflowequation(md,varargin) +% +% Example: +% md=setflowequation(md,'pattyn','Pattyn.exp','macayeal',md.mask.elementonfloatingice,'fill','hutter'); +% md=setflowequation(md,'pattyn','Pattyn.exp',fill','hutter','coupling','tiling'); + +%some checks on list of arguments +if ((nargin<2) | (nargout~=1)), + error('setflowequation error message'); +end + +%Find_out what kind of coupling to use +options=pairoptions(varargin{:}); +coupling_method=getfieldvalue(options,'coupling','tiling'); +if (~strcmpi(coupling_method,'tiling') & ~strcmpi(coupling_method,'penalties')), + error('coupling type can only be: tiling or penalties'); +end + +[hutterflag macayealflag pattynflag stokesflag filltype]=recover_areas(md,varargin{:}); + +%Flag the elements that have not been flagged as filltype +if strcmpi(filltype,'hutter'), + hutterflag(find(~(macayealflag | pattynflag)))=1; +elseif strcmpi(filltype,'macayeal'), + macayealflag(find(~(hutterflag | pattynflag | stokesflag)))=1; +elseif strcmpi(filltype,'pattyn'), + pattynflag(find(~(hutterflag | macayealflag | stokesflag)))=1; +end + +%check that each element has at least one flag +if any(hutterflag+ macayealflag+pattynflag+stokesflag==0), + error('setflowequation error message: elements type not assigned, must be specified') +end + +%check that each element has only one flag +if any(hutterflag+ macayealflag+pattynflag+stokesflag>1), + disp('setflowequation warning message: some elements have several types, higher order type is used for them') + hutterflag(find(hutterflag & macayealflag))=0; + hutterflag(find(hutterflag & pattynflag))=0; + macayealflag(find(macayealflag & pattynflag))=0; +end + +%Check that no pattyn or stokes for 2d mesh +if (md.mesh.dimension==2), + if any(stokesflag | pattynflag) + error('setflowequation error message: stokes and pattyn elements not allowed in 2d mesh, extrude it first') + end +end + +%Stokes can only be used alone for now: +if any(stokesflag) &any(hutterflag), + error('setflowequation error message: stokes cannot be used with any other model for now, put stokes everywhere') +end + +%Initialize node fields +nodeonhutter=zeros(md.mesh.numberofvertices,1); +nodeonhutter(md.mesh.elements(find(hutterflag),:))=1; +nodeonmacayeal=zeros(md.mesh.numberofvertices,1); +nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1; +nodeonpattyn=zeros(md.mesh.numberofvertices,1); +nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1; +nodeonstokes=zeros(md.mesh.numberofvertices,1); +noneflag=zeros(md.mesh.numberofelements,1); + +%First modify stokesflag to get rid of elements contrained everywhere (spc + border with pattyn or macayeal) +if any(stokesflag), + fullspcnodes=double((~isnan(md.diagnostic.spcvx)+~isnan(md.diagnostic.spcvy)+~isnan(md.diagnostic.spcvz))==3 | (nodeonpattyn & nodeonstokes)); %find all the nodes on the boundary of the domain without icefront + fullspcelems=double(sum(fullspcnodes(md.mesh.elements),2)==6); %find all the nodes on the boundary of the domain without icefront + stokesflag(find(fullspcelems))=0; + nodeonstokes(md.mesh.elements(find(stokesflag),:))=1; +end + +%Then complete with NoneApproximation or the other model used if there is no stokes +if any(stokesflag), + if any(pattynflag), %fill with pattyn + pattynflag(~stokesflag)=1; + nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1; + elseif any(macayealflag), %fill with macayeal + macayealflag(~stokesflag)=1; + nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1; + else %fill with none + noneflag(find(~stokesflag))=1; + end +end + +%Now take care of the coupling between MacAyeal and Pattyn +md.diagnostic.vertex_pairing=[]; +nodeonmacayealpattyn=zeros(md.mesh.numberofvertices,1); +nodeonpattynstokes=zeros(md.mesh.numberofvertices,1); +nodeonmacayealstokes=zeros(md.mesh.numberofvertices,1); +macayealpattynflag=zeros(md.mesh.numberofelements,1); +macayealstokesflag=zeros(md.mesh.numberofelements,1); +pattynstokesflag=zeros(md.mesh.numberofelements,1); +if strcmpi(coupling_method,'penalties'), + %Create the border nodes between Pattyn and MacAyeal and extrude them + numnodes2d=md.mesh.numberofvertices2d; + numlayers=md.mesh.numberoflayers; + bordernodes2d=find(nodeonpattyn(1:numnodes2d) & nodeonmacayeal(1:numnodes2d)); %Nodes connected to two different types of elements + + %initialize and fill in penalties structure + if ~isnan(bordernodes2d), + penalties=[]; + for i=1:numlayers-1, + penalties=[penalties; [bordernodes2d bordernodes2d+md.mesh.numberofvertices2d*(i)]]; + end + md.diagnostic.vertex_pairing=penalties; + end +elseif strcmpi(coupling_method,'tiling'), + if any(macayealflag) & any(pattynflag), %coupling macayeal pattyn + %Find node at the border + nodeonmacayealpattyn(find(nodeonmacayeal & nodeonpattyn))=1; + %Macayeal elements in contact with this layer become MacAyealPattyn elements + matrixelements=ismember(md.mesh.elements,find(nodeonmacayealpattyn)); + commonelements=sum(matrixelements,2)~=0; + commonelements(find(pattynflag))=0; %only one layer: the elements previously in macayeal + macayealflag(find(commonelements))=0; %these elements are now macayealpattynelements + macayealpattynflag(find(commonelements))=1; + nodeonmacayeal(:)=0; + nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1; + + %rule out elements that don't touch the 2 boundaries + pos=find(macayealpattynflag); + elist=zeros(length(pos),1); + elist = elist + any(sum(nodeonmacayeal(md.mesh.elements(pos,:)),2),2); + elist = elist - any(sum(nodeonpattyn(md.mesh.elements(pos,:)) ,2),2); + pos1=find(elist==1); + macayealflag(pos(pos1))=1; + macayealpattynflag(pos(pos1))=0; + pos2=find(elist==-1); + pattynflag(pos(pos2))=1; + macayealpattynflag(pos(pos2))=0; + + %Recompute nodes associated to these elements + nodeonmacayeal(:)=0; + nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1; + nodeonpattyn(:)=0; + nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1; + nodeonmacayealpattyn(:)=0; + nodeonmacayealpattyn(md.mesh.elements(find(macayealpattynflag),:))=1; + + elseif any(pattynflag) & any(stokesflag), %coupling pattyn stokes + %Find node at the border + nodeonpattynstokes(find(nodeonpattyn & nodeonstokes))=1; + %Stokes elements in contact with this layer become PattynStokes elements + matrixelements=ismember(md.mesh.elements,find(nodeonpattynstokes)); + commonelements=sum(matrixelements,2)~=0; + commonelements(find(pattynflag))=0; %only one layer: the elements previously in macayeal + stokesflag(find(commonelements))=0; %these elements are now macayealpattynelements + pattynstokesflag(find(commonelements))=1; + nodeonstokes=zeros(md.mesh.numberofvertices,1); + nodeonstokes(md.mesh.elements(find(stokesflag),:))=1; + + %rule out elements that don't touch the 2 boundaries + pos=find(pattynstokesflag); + elist=zeros(length(pos),1); + elist = elist + any(sum(nodeonstokes(md.mesh.elements(pos,:)),2),2); + elist = elist - any(sum(nodeonpattyn(md.mesh.elements(pos,:)),2),2); + pos1=find(elist==1); + stokesflag(pos(pos1))=1; + pattynstokesflag(pos(pos1))=0; + pos2=find(elist==-1); + pattynflag(pos(pos2))=1; + pattynstokesflag(pos(pos2))=0; + + %Recompute nodes associated to these elements + nodeonstokes(:)=0; + nodeonstokes(md.mesh.elements(find(stokesflag),:))=1; + nodeonpattyn(:)=0; + nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1; + nodeonpattynstokes(:)=0; + nodeonpattynstokes(md.mesh.elements(find(pattynstokesflag),:))=1; + + elseif any(stokesflag) & any(macayealflag), + %Find node at the border + nodeonmacayealstokes(find(nodeonmacayeal & nodeonstokes))=1; + %Stokes elements in contact with this layer become MacAyealStokes elements + matrixelements=ismember(md.mesh.elements,find(nodeonmacayealstokes)); + commonelements=sum(matrixelements,2)~=0; + commonelements(find(macayealflag))=0; %only one layer: the elements previously in macayeal + stokesflag(find(commonelements))=0; %these elements are now macayealmacayealelements + macayealstokesflag(find(commonelements))=1; + nodeonstokes=zeros(md.mesh.numberofvertices,1); + nodeonstokes(md.mesh.elements(find(stokesflag),:))=1; + + %rule out elements that don't touch the 2 boundaries + pos=find(macayealstokesflag); + elist=zeros(length(pos),1); + elist = elist + any(sum(nodeonmacayeal(md.mesh.elements(pos,:)),2),2); + elist = elist - any(sum(nodeonstokes(md.mesh.elements(pos,:)) ,2),2); + pos1=find(elist==1); + macayealflag(pos(pos1))=1; + macayealstokesflag(pos(pos1))=0; + pos2=find(elist==-1); + stokesflag(pos(pos2))=1; + macayealstokesflag(pos(pos2))=0; + + %Recompute nodes associated to these elements + nodeonmacayeal(:)=0; + nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1; + nodeonstokes(:)=0; + nodeonstokes(md.mesh.elements(find(stokesflag),:))=1; + nodeonmacayealstokes(:)=0; + nodeonmacayealstokes(md.mesh.elements(find(macayealstokesflag),:))=1; + + elseif any(stokesflag) & any(hutterflag), + error('type of coupling not supported yet'); + end +end + +%Create MacaAyealPattynApproximation where needed +md.flowequation.element_equation=zeros(md.mesh.numberofelements,1); +md.flowequation.element_equation(find(noneflag))=0; +md.flowequation.element_equation(find(hutterflag))=1; +md.flowequation.element_equation(find(macayealflag))=2; +md.flowequation.element_equation(find(pattynflag))=3; +md.flowequation.element_equation(find(stokesflag))=4; +md.flowequation.element_equation(find(macayealpattynflag))=5; +md.flowequation.element_equation(find(macayealstokesflag))=6; +md.flowequation.element_equation(find(pattynstokesflag))=7; + +%border +md.flowequation.borderpattyn=nodeonpattyn; +md.flowequation.bordermacayeal=nodeonmacayeal; +md.flowequation.borderstokes=nodeonstokes; + +%Create vertices_type +md.flowequation.vertex_equation=zeros(md.mesh.numberofvertices,1); +pos=find(nodeonhutter); +md.flowequation.vertex_equation(pos)=1; +pos=find(nodeonmacayeal); +md.flowequation.vertex_equation(pos)=2; +pos=find(nodeonpattyn); +md.flowequation.vertex_equation(pos)=3; +pos=find(nodeonhutter); +md.flowequation.vertex_equation(pos)=1; +pos=find(nodeonmacayealpattyn); +md.flowequation.vertex_equation(pos)=5; +pos=find(nodeonstokes); +md.flowequation.vertex_equation(pos)=4; +if any(stokesflag), + pos=find(~nodeonstokes); + if(~any(pattynflag) & ~any(macayealflag)), + md.flowequation.vertex_equation(pos)=0; + end +end +pos=find(nodeonpattynstokes); +md.flowequation.vertex_equation(pos)=7; +pos=find(nodeonmacayealstokes); +md.flowequation.vertex_equation(pos)=6; + +%figure out solution types +md.flowequation.ishutter=double(any(md.flowequation.element_equation==1)); +md.flowequation.ismacayealpattyn=double(any(md.flowequation.element_equation==2 | md.flowequation.element_equation==3)); +md.flowequation.isstokes=double(any(md.flowequation.element_equation==4)); + +return + +%Check that tiling can work: +if any(md.flowequation.bordermacayeal) & any(md.flowequation.borderpattyn) & any(md.flowequation.borderpattyn + md.flowequation.bordermacayeal ~=1), + error('error coupling domain too irregular'); +end +if any(md.flowequation.bordermacayeal) & any(md.flowequation.borderstokes) & any(md.flowequation.borderstokes + md.flowequation.bordermacayeal ~=1), + error('error coupling domain too irregular'); +end +if any(md.flowequation.borderstokes) & any(md.flowequation.borderpattyn) & any(md.flowequation.borderpattyn + md.flowequation.borderstokes~=1), + error('error coupling domain too irregular'); +end Index: /issm/trunk-jpl/src/m/model/parameterization/setflowequation.py =================================================================== --- /issm/trunk-jpl/src/m/model/parameterization/setflowequation.py (revision 13006) +++ /issm/trunk-jpl/src/m/model/parameterization/setflowequation.py (revision 13006) @@ -0,0 +1,278 @@ +import numpy +from model import * +from pairoptions import * +from recover_areas import * +from MatlabFuncs import * + +def setflowequation(md,*args): + """ + SETELEMENTSTYPE - associate a solution type to each element + + This routine works like plotmodel: it works with an even number of inputs + 'hutter','macayeal','pattyn','stokes' and 'fill' are the possible options + that must be followed by the corresponding exp file or flags list + It can either be a domain file (argus type, .exp extension), or an array of element flags. + If user wants every element outside the domain to be + setflowequationd, add '~' to the name of the domain file (ex: '~Pattyn.exp'); + an empty string '' will be considered as an empty domain + a string 'all' will be considered as the entire domain + You can specify the type of coupling, 'penalties' or 'tiling', to use with the input 'coupling' + + Usage: + md=setflowequation(md,varargin) + + Example: + md=setflowequation(md,'pattyn','Pattyn.exp','macayeal',md.mask.elementonfloatingice,'fill','hutter'); + md=setflowequation(md,'pattyn','Pattyn.exp',fill','hutter','coupling','tiling'); + """ + + #some checks on list of arguments + if not isinstance(md,model) or not len(args): + raise TypeError("setflowequation error message") + + #Find_out what kind of coupling to use + options=pairoptions(*args) + coupling_method=options.getfieldvalue('coupling','tiling') + if not strcmpi(coupling_method,'tiling') and not strcmpi(coupling_method,'penalties'): + raise TypeError("coupling type can only be: tiling or penalties") + + hutterflag,macayealflag,pattynflag,stokesflag,filltype=recover_areas(md,*args) + + #Flag the elements that have not been flagged as filltype + if strcmpi(filltype,'hutter'): + hutterflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(macayealflag,pattynflag)))]=1 + elif strcmpi(filltype,'macayeal'): + macayealflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(pattynflag,stokesflag))))]=1 + elif strcmpi(filltype,'pattyn'): + pattynflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(macayealflag,stokesflag))))]=1 + + #check that each element has at least one flag + if not numpy.any(hutterflag+macayealflag+pattynflag+stokesflag): + raise TypeError("setflowequation error message: elements type not assigned, must be specified") + + #check that each element has only one flag + if numpy.any(hutterflag+macayealflag+pattynflag+stokesflag>1): + print "setflowequation warning message: some elements have several types, higher order type is used for them" + hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,macayealflag))]=0 + hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,pattynflag))]=0 + macayealflag[numpy.nonzero(numpy.logical_and(macayealflag,pattynflag))]=0 + + #Check that no pattyn or stokes for 2d mesh + if md.mesh.dimension==2: + if numpy.any(numpy.logical_or(stokesflag,pattynflag)): + raise TypeError("setflowequation error message: stokes and pattyn elements not allowed in 2d mesh, extrude it first") + + #Stokes can only be used alone for now: + if numpy.any(stokesflag) and numpy.any(hutterflag): + raise TypeError("setflowequation error message: stokes cannot be used with any other model for now, put stokes everywhere") + + #Initialize node fields + nodeonhutter=numpy.zeros(md.mesh.numberofvertices) + nodeonhutter[md.mesh.elements[numpy.nonzero(hutterflag),:].astype(int)-1]=1 + nodeonmacayeal=numpy.zeros(md.mesh.numberofvertices) + nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1 + nodeonpattyn=numpy.zeros(md.mesh.numberofvertices) + nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1 + nodeonstokes=numpy.zeros(md.mesh.numberofvertices) + noneflag=numpy.zeros(md.mesh.numberofelements) + + #First modify stokesflag to get rid of elements contrained everywhere (spc + border with pattyn or macayeal) + if any(stokesflag): +# fullspcnodes=double((~isnan(md.diagnostic.spcvx)+~isnan(md.diagnostic.spcvy)+~isnan(md.diagnostic.spcvz))==3 | (nodeonpattyn & nodeonstokes)); %find all the nodes on the boundary of the domain without icefront + fullspcnodes=numpy.logical_or(numpy.logical_not(numpy.isnan(md.diagnostic.spcvx))+ \ + numpy.logical_not(numpy.isnan(md.diagnostic.spcvy))+ \ + numpy.logical_not(numpy.isnan(md.diagnostic.spcvz))==3, \ + numpy.logical_and(nodeonpattyn,nodeonstokes)).astype(int) #find all the nodes on the boundary of the domain without icefront +# fullspcelems=double(sum(fullspcnodes(md.mesh.elements),2)==6); %find all the nodes on the boundary of the domain without icefront + fullspcelems=(numpy.sum(fullspcnodes[md.mesh.elements.astype(int)-1],axis=1)==6).astype(int) #find all the nodes on the boundary of the domain without icefront + stokesflag[numpy.nonzero(fullspcelems)]=0 + nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1 + + #Then complete with NoneApproximation or the other model used if there is no stokes + if any(stokesflag): + if any(pattynflag): #fill with pattyn + pattynflag[numpy.logical_not(stokesflag)]=1 + nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1 + elif any(macayealflag): #fill with macayeal + macayealflag[numpy.logical_not(stokesflag)]=1 + nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1 + else: #fill with none + noneflag[numpy.nonzero(numpy.logical_not(stokesflag))]=1 + + #Now take care of the coupling between MacAyeal and Pattyn + md.diagnostic.vertex_pairing=numpy.array([]) + nodeonmacayealpattyn=numpy.zeros(md.mesh.numberofvertices) + nodeonpattynstokes=numpy.zeros(md.mesh.numberofvertices) + nodeonmacayealstokes=numpy.zeros(md.mesh.numberofvertices) + macayealpattynflag=numpy.zeros(md.mesh.numberofelements) + macayealstokesflag=numpy.zeros(md.mesh.numberofelements) + pattynstokesflag=numpy.zeros(md.mesh.numberofelements) + if strcmpi(coupling_method,'penalties'): + #Create the border nodes between Pattyn and MacAyeal and extrude them + numnodes2d=md.mesh.numberofvertices2d + numlayers=md.mesh.numberoflayers + bordernodes2d=numpy.nonzero(numpy.logical_and(nodeonpattyn[1:numnodes2d],nodeonmacayeal[1:numnodes2d])) #Nodes connected to two different types of elements + + #initialize and fill in penalties structure + if numpy.all(numpy.logical_not(numpy.isnan(bordernodes2d))): + penalties=numpy.zeros((0,2)) + for i in xrange(1,numlayers): + penalties=numpy.concatenate((penalties,numpy.concatenate((bordernodes2d,bordernodes2d+md.mesh.numberofvertices2d*(i)),axis=1)),axis=0) + md.diagnostic.vertex_pairing=penalties + + elif strcmpi(coupling_method,'tiling'): + if numpy.any(macayealflag) and numpy.any(pattynflag): #coupling macayeal pattyn + #Find node at the border + nodeonmacayealpattyn[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonpattyn))]=1 + #Macayeal elements in contact with this layer become MacAyealPattyn elements + matrixelements=ismember(md.mesh.elements,numpy.nonzero(nodeonmacayealpattyn)) + commonelements=numpy.sum(matrixelements,axis=1)!=0 + commonelements[numpy.nonzero(pattynflag)]=0 #only one layer: the elements previously in macayeal + macayealflag[numpy.nonzero(commonelements)]=0 #these elements are now macayealpattynelements + macayealpattynflag[numpy.nonzero(commonelements)]=1 + nodeonmacayeal[:]=0 + nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1 + + #rule out elements that don't touch the 2 boundaries + pos=numpy.nonzero(macayealpattynflag) + elist=numpy.zeros(len(pos)) + elist = elist + numpy.any(numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:].astype(int)-1],axis=1),axis=1) + elist = elist - numpy.any(numpy.sum(nodeonpattyn[md.mesh.elements[pos,:].astype(int)-1] ,axis=1),axis=1) + pos1=[i for i,item in enumerate(elist) if item==1] + macayealflag[pos[pos1]]=1 + macayealpattynflag[pos[pos1]]=0 + pos2=[i for i,item in enumerate(elist) if item==-1] + pattynflag[pos[pos2]]=1 + macayealpattynflag[pos[pos2]]=0 + + #Recompute nodes associated to these elements + nodeonmacayeal[:]=0 + nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1 + nodeonpattyn[:]=0 + nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1 + nodeonmacayealpattyn[:]=0 + nodeonmacayealpattyn[md.mesh.elements[numpy.nonzero(macayealpattynflag),:].astype(int)-1]=1 + + elif numpy.any(pattynflag) and numpy.any(stokesflag): #coupling pattyn stokes + #Find node at the border + nodeonpattynstokes[numpy.nonzero(numpy.logical_and(nodeonpattyn,nodeonstokes))]=1 + #Stokes elements in contact with this layer become PattynStokes elements + matrixelements=ismember(md.mesh.elements,numpy.nonzero(nodeonpattynstokes)) + commonelements=numpy.sum(matrixelements,axis=1)!=0 + commonelements[numpy.nonzero(pattynflag)]=0 #only one layer: the elements previously in macayeal + stokesflag[numpy.nonzero(commonelements)]=0 #these elements are now macayealpattynelements + pattynstokesflag[numpy.nonzero(commonelements)]=1 + nodeonstokes=numpy.zeros(md.mesh.numberofvertices) + nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1 + + #rule out elements that don't touch the 2 boundaries + pos=numpy.nonzero(pattynstokesflag) + elist=numpy.zeros(len(pos)) + elist = elist + numpy.any(numpy.sum(nodeonstokes[md.mesh.elements[pos,:].astype(int)-1],axis=1),axis=1) + elist = elist - numpy.any(numpy.sum(nodeonpattyn[md.mesh.elements[pos,:].astype(int)-1],axis=1),axis=1) + pos1=[i for i,item in enumerate(elist) if item==1] + stokesflag[pos[pos1]]=1 + pattynstokesflag[pos[pos1]]=0 + pos2=[i for i,item in enumerate(elist) if item==-1] + pattynflag[pos[pos2]]=1 + pattynstokesflag[pos[pos2]]=0 + + #Recompute nodes associated to these elements + nodeonstokes[:]=0 + nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1 + nodeonpattyn[:]=0 + nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1 + nodeonpattynstokes[:]=0 + nodeonpattynstokes[md.mesh.elements[numpy.nonzero(pattynstokesflag),:].astype(int)-1]=1 + + elif numpy.any(stokesflag) and numpy.any(macayealflag): + #Find node at the border + nodeonmacayealstokes[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonstokes))]=1 + #Stokes elements in contact with this layer become MacAyealStokes elements + matrixelements=ismember(md.mesh.elements,numpy.nonzero(nodeonmacayealstokes)) + commonelements=numpy.sum(matrixelements,axis=1)!=0 + commonelements[numpy.nonzero(macayealflag)]=0 #only one layer: the elements previously in macayeal + stokesflag[numpy.nonzero(commonelements)]=0 #these elements are now macayealmacayealelements + macayealstokesflag[numpy.nonzero(commonelements)]=1 + nodeonstokes=numpy.zeros(md.mesh.numberofvertices) + nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1 + + #rule out elements that don't touch the 2 boundaries + pos=numpy.nonzero(macayealstokesflag) + elist=numpy.zeros(len(pos)) + elist = elist + numpy.any(numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:].astype(int)-1],axis=1),axis=1) + elist = elist - numpy.any(numpy.sum(nodeonstokes[md.mesh.elements[pos,:].astype(int)-1] ,axis=1),axis=1) + pos1=[i for i,item in enumerate(elist) if item==1] + macayealflag[pos[pos1]]=1 + macayealstokesflag[pos[pos1]]=0 + pos2=[i for i,item in enumerate(elist) if item==-1] + stokesflag[pos[pos2]]=1 + macayealstokesflag[pos[pos2]]=0 + + #Recompute nodes associated to these elements + nodeonmacayeal[:]=0 + nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1 + nodeonstokes[:]=0 + nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1 + nodeonmacayealstokes[:]=0 + nodeonmacayealstokes[md.mesh.elements[numpy.nonzero(macayealstokesflag),:].astype(int)-1]=1 + + elif numpy.any(stokesflag) and numpy.any(hutterflag): + raise TypeError("type of coupling not supported yet") + + #Create MacaAyealPattynApproximation where needed + md.flowequation.element_equation=numpy.zeros(md.mesh.numberofelements) + md.flowequation.element_equation[numpy.nonzero(noneflag)]=0 + md.flowequation.element_equation[numpy.nonzero(hutterflag)]=1 + md.flowequation.element_equation[numpy.nonzero(macayealflag)]=2 + md.flowequation.element_equation[numpy.nonzero(pattynflag)]=3 + md.flowequation.element_equation[numpy.nonzero(stokesflag)]=4 + md.flowequation.element_equation[numpy.nonzero(macayealpattynflag)]=5 + md.flowequation.element_equation[numpy.nonzero(macayealstokesflag)]=6 + md.flowequation.element_equation[numpy.nonzero(pattynstokesflag)]=7 + + #border + md.flowequation.borderpattyn=nodeonpattyn + md.flowequation.bordermacayeal=nodeonmacayeal + md.flowequation.borderstokes=nodeonstokes + + #Create vertices_type + md.flowequation.vertex_equation=numpy.zeros(md.mesh.numberofvertices) + pos=numpy.nonzero(nodeonhutter) + md.flowequation.vertex_equation[pos]=1 + pos=numpy.nonzero(nodeonmacayeal) + md.flowequation.vertex_equation[pos]=2 + pos=numpy.nonzero(nodeonpattyn) + md.flowequation.vertex_equation[pos]=3 + pos=numpy.nonzero(nodeonhutter) + md.flowequation.vertex_equation[pos]=1 + pos=numpy.nonzero(nodeonmacayealpattyn) + md.flowequation.vertex_equation[pos]=5 + pos=numpy.nonzero(nodeonstokes) + md.flowequation.vertex_equation[pos]=4 + if numpy.any(stokesflag): + pos=numpy.nonzero(numpy.logical_not(nodeonstokes)) + if not (numpy.any(pattynflag) or numpy.any(macayealflag)): + md.flowequation.vertex_equation[pos]=0 + pos=numpy.nonzero(nodeonpattynstokes) + md.flowequation.vertex_equation[pos]=7 + pos=numpy.nonzero(nodeonmacayealstokes) + md.flowequation.vertex_equation[pos]=6 + + #figure out solution types + md.flowequation.ishutter=float(numpy.any(md.flowequation.element_equation==1)) + md.flowequation.ismacayealpattyn=float(numpy.any(numpy.logical_or(md.flowequation.element_equation==2,md.flowequation.element_equation==3))) + md.flowequation.isstokes=float(numpy.any(md.flowequation.element_equation==4)) + + return md + + #Check that tiling can work: + if numpy.any(md.flowequation.bordermacayeal) and numpy.any(md.flowequation.borderpattyn) and numpy.any(md.flowequation.borderpattyn + md.flowequation.bordermacayeal !=1): + raise TypeError("error coupling domain too irregular") + if numpy.any(md.flowequation.bordermacayeal) and numpy.any(md.flowequation.borderstokes) and numpy.any(md.flowequation.borderstokes + md.flowequation.bordermacayeal !=1): + raise TypeError("error coupling domain too irregular") + if numpy.any(md.flowequation.borderstokes) and numpy.any(md.flowequation.borderpattyn) and numpy.any(md.flowequation.borderpattyn + md.flowequation.borderstokes !=1): + raise TypeError("error coupling domain too irregular") + + return md + Index: /issm/trunk-jpl/src/m/model/parameterization/setmask.m =================================================================== --- /issm/trunk-jpl/src/m/model/parameterization/setmask.m (revision 13006) +++ /issm/trunk-jpl/src/m/model/parameterization/setmask.m (revision 13006) @@ -0,0 +1,49 @@ +function md=setmask(md,floatingicename,groundedicename) +%SETMASK - establish boundaries between grounded and floating ice. +% +% By default, ice is considered grounded. The contour floatingicename defines nodes +% for which ice is floating. The contour groundedicename defines nodes inside an floatingice, +% that are grounded (ie: ice rises, islands, etc ...) +% All input files are in the Argus format (extension .exp). +% +% Usage: +% md=setmask(md,floatingicename,groundedicename) +% +% Examples: +% md=setmask(md,'all',''); +% md=setmask(md,'Iceshelves.exp','Islands.exp'); + +%some checks on list of arguments +if ((nargin~=3) | (nargout~=1)), + help mask + error('mask error message'); +end + +%Get assigned fields +x=md.mesh.x; +y=md.mesh.y; +elements=md.mesh.elements; + +%Assign elementonfloatingice, elementongroundedice, vertexongroundedice and vertexonfloatingice. Only change at your own peril! This is synchronized heavily with the GroundingLineMigration module. {{{ +elementonfloatingice=FlagElements(md,floatingicename); +elementongroundedice=FlagElements(md,groundedicename); + +%Because groundedice nodes and elements can be included into an floatingice, we need to update. Remember, all the previous +%arrays come from domain outlines that can intersect one another: +elementonfloatingice=double((elementonfloatingice & ~elementongroundedice)); +elementongroundedice=double(~elementonfloatingice); + +%the order here is important. we choose vertexongroundedice as default on the grounding line. +vertexonfloatingice=zeros(md.mesh.numberofvertices,1); +vertexongroundedice=zeros(md.mesh.numberofvertices,1); +vertexongroundedice(md.mesh.elements(find(elementongroundedice),:))=1; +vertexonfloatingice(find(~vertexongroundedice))=1; +%}}} + +%Return: +md.mask.elementonfloatingice=elementonfloatingice; +md.mask.vertexonfloatingice=vertexonfloatingice; +md.mask.elementongroundedice=elementongroundedice; +md.mask.vertexongroundedice=vertexongroundedice; +md.mask.vertexonwater=zeros(md.mesh.numberofvertices,1); +md.mask.elementonwater=zeros(md.mesh.numberofelements,1); Index: /issm/trunk-jpl/src/m/model/parameterization/setmask.py =================================================================== --- /issm/trunk-jpl/src/m/model/parameterization/setmask.py (revision 13006) +++ /issm/trunk-jpl/src/m/model/parameterization/setmask.py (revision 13006) @@ -0,0 +1,55 @@ +from numpy import * +import FlagElements as fe + +def setmask(md, floatingicename, groundedicename): + #SETMASK - establish boundaries between grounded and floating ice. + # + # By default, ice is considered grounded. The contour floatingicename defines nodes + # for which ice is floating. The contour groundedicename defines nodes inside an floatingice, + # that are grounded (ie: ice rises, islands, etc ...) + # All input files are in the Argus format (extension .exp). + # + # Usage: + # md=setmask(md,floatingicename,groundedicename) + # + # Examples: + # md=setmask(md,'all',''); + # md=setmask(md,'Iceshelves.exp','Islands.exp'); + + #%Get assigned fields + x = md.mesh.x + y = md.mesh.y + elements = md.mesh.elements + + #Assign elementonfloatingice, elementongroundedice, vertexongroundedice and vertexonfloatingice. Only change at your own peril! This is synchronized heavily with the GroundingLineMigration module. {{{ + elementonfloatingice = fe.FlagElements(md, floatingicename) + elementongroundedice = fe.FlagElements(md, groundedicename) + + #Because groundedice nodes and elements can be included into an floatingice, we need to update. Remember, all the previous + #arrays come from domain outlines that can intersect one another: + + elementonfloatingice = logical_and(elementonfloatingice,~elementongroundedice) + elementongroundedice = ~elementonfloatingice + + #the order here is important. we choose vertexongroundedice as default on the grounding line. + vertexonfloatingice = zeros(md.mesh.numberofvertices,'bool') + vertexongroundedice = zeros(md.mesh.numberofvertices,'bool') + + pos=argwhere(elementongroundedice==1) + pos=md.mesh.elements[pos,:]-1 + if pos.size: + vertexongroundedice[pos]=True + + pos=argwhere(~vertexongroundedice) + if pos.size: + vertexonfloatingice[pos]=True; + #%}}} + + #Return: + md.mask.elementonfloatingice = double(elementonfloatingice) + md.mask.vertexonfloatingice = double(vertexonfloatingice) + md.mask.elementongroundedice = double(elementongroundedice) + md.mask.vertexongroundedice = double(vertexongroundedice) + md.mask.vertexonwater = zeros(md.mesh.numberofvertices) + md.mask.elementonwater = zeros(md.mesh.numberofelements) + return md Index: /issm/trunk-jpl/src/m/model/parameterization/setmask2.m =================================================================== --- /issm/trunk-jpl/src/m/model/parameterization/setmask2.m (revision 13006) +++ /issm/trunk-jpl/src/m/model/parameterization/setmask2.m (revision 13006) @@ -0,0 +1,148 @@ +function md=setmask2(md,landname,floatingicename,groundedicename) +%GEOGRAPHY2 - establish land, ice sheet and ice shelf areas in a domains. +% +% Usage: +% md=setmask2(md,landname,floatingicename,groundedicename) +% +% Examples: +% md=setmask2(md,'LandName.exp','Iceshelves.exp','Islands.exp'); + +%Get assigned fields +x=md.mesh.x; +y=md.mesh.y; +elements=md.mesh.elements; + +%recover elements and nodes on land. +if ischar(landname), + [vertexonland,elementonland]=ContourToMesh(elements,x,y,landname,'element and node',2); +elseif isfloat(landname), + if size(landname,1)~=md.mesh.numberofelements, + error('Landname for area must be of same size as number of elements in model'); + end + elementonland=landname; + vertexonland=zeros(md.mesh.numberofvertices,1); + vertexonland(md.mesh.elements(find(elementonland),:))=1; +else + error('Invalid area option option'); +end + +%Now, build the connectivity tables for this mesh. +if size(md.mesh.vertexconnectivity,1)~=md.mesh.numberofvertices, + md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices); +end +if size(md.mesh.elementconnectivity,1)~=md.mesh.numberofelements, + md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity); +end + +%any element with 3 nodes on land should be on land: +elementsonwater=find(~elementonland); +wrongelements=elementsonwater(find(( vertexonland(md.mesh.elements(elementsonwater,1)) + vertexonland(md.mesh.elements(elementsonwater,2)) + vertexonland(md.mesh.elements(elementsonwater,3)) ... + )==3)); +elementonland(wrongelements)=1; + +%any element with its barycentre on land should be on land: (only if landname is an expfile) +if ischar(landname), +weights={[1;1;1],[2;1;1],[1;2;1],[1;1;2]}; + for i=1:length(weights), + xelem=x(md.mesh.elements)*weights{i}/sum(weights{i}); + yelem=y(md.mesh.elements)*weights{i}/sum(weights{i}); + end + baryonland=ContourToNodes(xelem,yelem,landname,1); + pos=find(~baryonland); elementonland(pos)=0; + pos=find(baryonland); elementonland(pos)=1; +end + +%figure out which elements on land are actually in the middle of the ocean! +pos1=find(elementonland); +connectedtoland=md.mesh.elementconnectivity(pos1,:); +pos=find(connectedtoland); connectedtoland(pos)=1-elementonland(connectedtoland(pos)); +connectedtolandsum=sum(connectedtoland,2); +waterelements=pos1(find(connectedtolandsum==3)); +elementonland(waterelements)=0; + +%figure out which elements on water are actually in the middle of the land! +pos1=find(~elementonland); +connectedtowater=md.mesh.elementconnectivity(pos1,:); +pos=find(connectedtowater); connectedtowater(pos)=elementonland(connectedtowater(pos)); +connectedtowatersum=sum(connectedtowater,2); +landelements=pos1(find(connectedtowatersum==3)); +elementonland(landelements)=1; + +%recover arrays of ice shelf nodes and elements, and ice sheet nodes and elements. +elementonfloatingice=FlagElements(md,floatingicename); +elementongroundedice=FlagElements(md,groundedicename); + +%Because groundedice nodes and elements can be included into an floatingice, we need to update. Remember, all the previous +%arrays come from domain outlines that can intersect one another: +vertexonfloatingice=zeros(md.mesh.numberofvertices,1); +vertexongroundedice=zeros(md.mesh.numberofvertices,1); +elementonfloatingice=double((elementonfloatingice & ~elementongroundedice)); +elementongroundedice=double(~elementonfloatingice); +vertexonfloatingice(md.mesh.elements(find(elementonfloatingice),:))=1; +vertexongroundedice(md.mesh.elements(find(elementongroundedice),:))=1; + +%now correct, so that none of the floatingice and groundedice elements and nodes are in the water. +pos=find(~elementonland); +elementonfloatingice(pos)=0; +elementongroundedice(pos)=0; + +pos=find(~vertexonland); +vertexonfloatingice(pos)=0; +vertexongroundedice(pos)=0; + +%create vertexonwater and elementonwater: +vertexonwater=double(~vertexonland); +elementonwater=double(~elementonland); + +%correct for islands: +vertexonfloatingice=double(vertexonfloatingice & ~vertexongroundedice); +elementonfloatingice=double(elementonfloatingice & ~elementongroundedice); + +%now, groundedices are everything except iceshelves and water +vertexongroundedice=double(~vertexonfloatingice & ~vertexonwater); +elementongroundedice=double(~elementonfloatingice & ~elementonwater); + +%Deal with segments on neumann: + +%Get current connectivity +mesh.elementconnectivity=md.mesh.elementconnectivity; + +%put 0 for elements on water +pos=find(mesh.elementconnectivity); +mesh.elementconnectivity(pos)=mesh.elementconnectivity(pos).*(~elementonwater(mesh.elementconnectivity(pos))); + +%put line of ones for elements on water +pos=find(elementonwater); +mesh.elementconnectivity(pos,:)=1;% line of ones for elements on water so they won't be considered + +%resort lines (zeros must be at the last column for findsegments) +mesh.elementconnectivity=sort(mesh.elementconnectivity,2,'descend'); + +%call findsegments to build segment using THIS conectivity +md.mesh.segments=findsegments(md,'mesh.elementconnectivity',mesh.elementconnectivity); + +%some final checks: +%check that no node thinks it's on an ice shelf or ice sheet, and lies actually in the middle of the water. +nodesgrounded=find(~vertexonwater); +lengthconnectivity=size(md.mesh.vertexconnectivity,2); +groundedcounters=md.mesh.vertexconnectivity(nodesgrounded,lengthconnectivity); +groundedconnectivity=md.mesh.vertexconnectivity(nodesgrounded,1:lengthconnectivity-1); +pos=find(groundedconnectivity); +groundedconnectivity(pos)=elementonwater(groundedconnectivity(pos)); +groundedsum=sum(groundedconnectivity,2); +errorflags=find(groundedsum==groundedcounters); +errornodes=nodesgrounded(errorflags); + +vertexonwater(errornodes)=1; +vertexongroundedice(errornodes)=0; +vertexonfloatingice(errornodes)=0; + +%Return: +md.mask.vertexonfloatingice=vertexonfloatingice; +md.mask.elementonfloatingice=elementonfloatingice; +md.mask.vertexonwater=vertexonwater; +md.mask.elementonwater=elementonwater; +md.mask.vertexongroundedice=vertexongroundedice; +md.mask.elementongroundedice=elementongroundedice; + +md.mesh.segmentmarkers(:)=1; Index: sm/trunk-jpl/src/m/model/parseresultsfromdisk.m =================================================================== --- /issm/trunk-jpl/src/m/model/parseresultsfromdisk.m (revision 13005) +++ (revision ) @@ -1,212 +1,0 @@ -function results=parseresultsfromdisk(filename,iosplit) -%PARSERESULTSFROMDISK - ... -% -% Usage: -% results=parseresultsfromdisk(filename,iosplit) - -if iosplit, - results=parseresultsfromdiskiosplit(filename); -else - results=parseresultsfromdiskioserial(filename); -end - -%process patch if necessary -results=MatlabProcessPatch(results); - -function results=parseresultsfromdiskioserial(filename) % {{{ -%PARSERESULTSFROMDISK - ... -% -% Usage: -% results=parseresultsfromdiskioserial(filename) - - -%Open file -fid=fopen(filename,'rb'); -if(fid==-1), - error(['loadresultsfromdisk error message: could not open ',filename,' for binary reading']); -end -results=struct(); - -%Read fields until the end of the file. -result=ReadData(fid); -while ~isempty(result), - %Get time and step - results(result.step).step=result.step; - results(result.step).time=result.time; - - %Add result - if (length(results)>=result.step & isfield(results,result.fieldname) & ~strcmp(result.fieldname,'SolutionType')), - results(result.step).(result.fieldname)=[ results(result.step).(result.fieldname); result.field]; - else - results(result.step).(result.fieldname)=result.field; - end - - %read next result - result=ReadData(fid); - -end - -fclose(fid); -% }}} -function results=parseresultsfromdiskiosplit(filename) % {{{ -%PARSERESULTSFROMDISKIOSPLIT - ... -% -% Usage: -% results=parseresultsfromdiskiosplit(filename) - - -%Open file -fid=fopen(filename,'rb'); -if(fid==-1), - error(['loadresultsfromdisk error message: could not open ',filename,' for binary reading']); -end -results=struct(); - -%if we have done split I/O, ie, we have results that are fragmented across patches, -%do a first pass, and figure out the structure of results -result=ReadDataDimensions(fid); -while ~isempty(result), - - %Get time and step - results(result.step).step=result.step; - results(result.step).time=result.time; - - %Add result - if strcmpi(result.fieldname,'Patch'), - results(result.step).(result.fieldname)=[0 result.N]; - else - results(result.step).(result.fieldname)=NaN; - end - - %read next result - result=ReadDataDimensions(fid); -end - -%do a second pass, and figure out the size of the patches -fseek(fid,0,-1); %rewind -result=ReadDataDimensions(fid); -while ~isempty(result), - - %Add result - if strcmpi(result.fieldname,'Patch'), - patchdimensions=results(result.step).(result.fieldname); - results(result.step).(result.fieldname)=[patchdimensions(1)+result.M result.N]; - end - - %read next result - result=ReadDataDimensions(fid); -end - -%allocate patches -for i=1:length(results), - results(i).Patch=zeros(results(i).Patch(1),results(i).Patch(2)); - results(i).counter=1; %use to index into the patch -end - -%third pass, this time to read the real information -fseek(fid,0,-1); %rewind -result=ReadData(fid); -while ~isempty(result), - - %Get time and step - results(result.step).step=result.step; - results(result.step).time=result.time; - - %Add result - if strcmpi(result.fieldname,'Patch'), - counter=results(result.step).counter; - counter2=counter+size(result.field,1)-1; - results(result.step).(result.fieldname)(counter:counter2,:)=result.field; - - %increment counter: - results(result.step).counter=counter2+1; - else - results(result.step).(result.fieldname)=result.field; - end - - %read next result - result=ReadData(fid); - -end - -%close file -fclose(fid); - % }}} -function result=ReadData(fid) % {{{ -%READDATA - ... -% -% Usage: -% field=ReadData(fid) - -%read field -[length,count]=fread(fid,1,'int'); - -if count==0, - result=struct([]); -else - fieldname=fread(fid,length,'char'); - fieldname=fieldname(1:end-1)'; - fieldname=char(fieldname); - time=fread(fid,1,'double'); - step=fread(fid,1,'int'); - - type=fread(fid,1,'int'); - M=fread(fid,1,'int'); - if type==1, - field=fread(fid,M,'double'); - elseif type==2, - field=fread(fid,M,'char'); - field=char(field(1:end-1)'); - elseif type==3, - N=fread(fid,1,'int'); - field=transpose(fread(fid,[N M],'double')); - else - error(['cannot read data of type ' num2str(type) ]); - end - - result.fieldname=fieldname; - result.time=time; - result.step=step; - result.field=field; -end -% }}} -function result=ReadDataDimensions(fid) % {{{ -%READDATADIMENSIONS - read data dimensions, step and time, but not the data itself. -% -% Usage: -% field=ReadDataDimensions(fid) - - -%read field -[length,count]=fread(fid,1,'int'); - -if count==0, - result=struct([]); -else - fieldname=fread(fid,length,'char'); - fieldname=fieldname(1:end-1)'; - fieldname=char(fieldname); - time=fread(fid,1,'double'); - step=fread(fid,1,'int'); - - type=fread(fid,1,'int'); - M=fread(fid,1,'int'); - N=1; %default - if type==1, - fseek(fid,M*8,0); - elseif type==2, - fseek(fid,M,0); - elseif type==3, - N=fread(fid,1,'int'); - fseek(fid,N*M*8,0); - else - error(['cannot read data of type ' num2str(type) ]); - end - - result.fieldname=fieldname; - result.time=time; - result.step=step; - result.M=M; - result.N=N; -end -% }}} Index: sm/trunk-jpl/src/m/model/parseresultsfromdisk.py =================================================================== --- /issm/trunk-jpl/src/m/model/parseresultsfromdisk.py (revision 13005) +++ (revision ) @@ -1,235 +1,0 @@ -import struct -import numpy -from MatlabFuncs import * -from MatlabProcessPatch import * - -def parseresultsfromdisk(filename,iosplit): - """ - PARSERESULTSFROMDISK - ... - - Usage: - results=parseresultsfromdisk(filename,iosplit) - """ - - if iosplit: - results=parseresultsfromdiskiosplit(filename) - else: - results=parseresultsfromdiskioserial(filename) - - #process patch if necessary - results=MatlabProcessPatch(results) - - return results - -def parseresultsfromdiskioserial(filename): # {{{ - """ - PARSERESULTSFROMDISK - ... - - Usage: - results=parseresultsfromdiskioserial(filename) - """ - - #Open file - try: - fid=open(filename,'rb') - except IOError as e: - raise IOError("loadresultsfromdisk error message: could not open '%s' for binary reading." % filename) - - results={} - - #Read fields until the end of the file. - result=ReadData(fid) - while result: - #Get time and step - if not result['step'] in results: - results[result['step']]={} - results[result['step']]['step']=result['step'] - results[result['step']]['time']=result['time'] - - #Add result - if result['step'] in results and \ - result['fieldname'] in results[result['step']] and \ - not strcmp(result['fieldname'],'SolutionType'): - results[result['step']][result['fieldname']]=numpy.concatenate((results[result['step']][result['fieldname']],result['field']),axis=0) - else: - results[result['step']][result['fieldname']]=result['field'] - - #read next result - result=ReadData(fid) - - fid.close() - - return results - # }}} - -def parseresultsfromdiskiosplit(filename): # {{{ - """ - PARSERESULTSFROMDISKIOSPLIT - ... - - Usage: - results=parseresultsfromdiskiosplit(filename) - """ - - #Open file - try: - fid=open(filename,'rb') - except IOError as e: - raise IOError("loadresultsfromdisk error message: could not open '%s' for binary reading." % filename) - - results={} - - #if we have done split I/O, ie, we have results that are fragmented across patches, - #do a first pass, and figure out the structure of results - result=ReadDataDimensions(fid) - while result: - - #Get time and step - if not result['step'] in results: - results[result['step']]={} - results[result['step']]['step']=result['step'] - results[result['step']]['time']=result['time'] - - #Add result - if strcmpi(result['fieldname'],'Patch'): - results[result['step']][result['fieldname']]=[0,result['N']] - else: - results[result['step']][result['fieldname']]=float('NaN') - - #read next result - result=ReadDataDimensions(fid) - - #do a second pass, and figure out the size of the patches - fid.seek(0) #rewind - result=ReadDataDimensions(fid) - while result: - - #Add result - if strcmpi(result['fieldname'],'Patch'): - patchdimensions=results[result['step']][result['fieldname']] - results[result['step']][result['fieldname']]=[patchdimensions[0]+result['M'],result['N']] - - #read next result - result=ReadDataDimensions(fid) - - #allocate patches - for result in results.itervalues(): - if 'Patch' in result: - result['Patch']=numpy.zeros(shape=(result['Patch'][0],result['Patch'][1]),dtype=float) - result['counter']=0 #use to index into the patch - - #third pass, this time to read the real information - fid.seek(0) #rewind - result=ReadData(fid) - while result: - - #Get time and step - if not result['step'] in results: - results[result['step']]={} - results[result['step']]['step']=result['step'] - results[result['step']]['time']=result['time'] - - #Add result - if strcmpi(result['fieldname'],'Patch'): - counter=results[result['step']]['counter'] - counter2=counter+result['field'].shape[0]-1 - results[result['step']][result['fieldname']][counter:counter2,:]=result['field'] - - #increment counter: - results[result['step']]['counter']=counter2+1 - else: - results[result['step']][result['fieldname']]=result['field'] - - #read next result - result=ReadData(fid) - - #close file - fid.close() - - return results - # }}} - -def ReadData(fid): # {{{ - """ - READDATA - ... - - Usage: - field=ReadData(fid) - """ - - #read field - try: - length=struct.unpack('i',fid.read(struct.calcsize('i')))[0] - - fieldname=struct.unpack('%ds' % length,fid.read(length))[0][:-1] - time=struct.unpack('d',fid.read(struct.calcsize('d')))[0] - step=struct.unpack('i',fid.read(struct.calcsize('i')))[0] - - type=struct.unpack('i',fid.read(struct.calcsize('i')))[0] - M=struct.unpack('i',fid.read(struct.calcsize('i')))[0] - if type==1: - field=numpy.array(struct.unpack('%dd' % M,fid.read(M*struct.calcsize('d'))),dtype=float) - elif type==2: - field=struct.unpack('%ds' % M,fid.read(M))[0][:-1] - elif type==3: - N=struct.unpack('i',fid.read(struct.calcsize('i')))[0] -# field=transpose(fread(fid,[N M],'double')); - field=numpy.zeros(shape=(M,N),dtype=float) - for i in xrange(M): - field[i,:]=struct.unpack('%dd' % N,fid.read(N*struct.calcsize('d'))) - else: - raise TypeError("cannot read data of type %d" % type) - - result={} - result['fieldname']=fieldname - result['time']=time - result['step']=step - result['field']=field - - except struct.error as e: - result={} - - return result - # }}} - -def ReadDataDimensions(fid): # {{{ - """ - READDATADIMENSIONS - read data dimensions, step and time, but not the data itself. - - Usage: - field=ReadDataDimensions(fid) - """ - - #read field - try: - length=struct.unpack('i',fid.read(struct.calcsize('i')))[0] - - fieldname=struct.unpack('%ds' % length,fid.read(length))[0][:-1] - time=struct.unpack('d',fid.read(struct.calcsize('d')))[0] - step=struct.unpack('i',fid.read(struct.calcsize('i')))[0] - - type=struct.unpack('i',fid.read(struct.calcsize('i')))[0] - M=struct.unpack('i',fid.read(struct.calcsize('i')))[0] - N=1 #default - if type==1: - fid.seek(M*8,1) - elif type==2: - fid.seek(M,1) - elif type==3: - N=struct.unpack('i',fid.read(struct.calcsize('i')))[0] - fid.seek(N*M*8,1) - else: - raise TypeError("cannot read data of type %d" % type) - - result={} - result['fieldname']=fieldname - result['time']=time - result['step']=step - result['M']=M - result['N']=N - - except struct.error as e: - result={} - - return result - # }}} - Index: /issm/trunk-jpl/src/m/model/print/printmodel.m =================================================================== --- /issm/trunk-jpl/src/m/model/print/printmodel.m (revision 13006) +++ /issm/trunk-jpl/src/m/model/print/printmodel.m (revision 13006) @@ -0,0 +1,100 @@ +function printmodel(filename,format,varargin) +%PRINTMODEL - save an image of current figure +% +% filename: output name of image file (no extension) +% format: image format (ex: 'tiff','jpg','pdf') +% +% List of options to printfmodel: +% +% figure: number of figure to print (default: current figure) +% resolution: use higher resolution to anti-alias (default 2) +% margin: add margin around final image +% marginsize: size of margin around final image (default 5) +% frame: add frame around final image +% framesize: size of frame around final image (default 5) +% framecolor: color of frame around final image (default 'black') +% trim: trim empty space around image (default 'off') +% hardcopy: 'off' to impose MATLAB to use the same colors (default 'off') +% +% Usage: +% printmodel(filename,format,varargin); +% +% Examples: +% printmodel('image','tiff') +% printmodel('image','eps','margin','on','frame','on','hardcopy','on') + + +%get options: +options=pairoptions(varargin{:}); + +%set defaults +options=addfielddefault(options,'figure','gcf'); +options=addfielddefault(options,'format','tiff'); +options=addfielddefault(options,'resolution',1); +options=addfielddefault(options,'margin','on'); +options=addfielddefault(options,'marginsize',25); +options=addfielddefault(options,'frame','on'); +options=addfielddefault(options,'framesize',3); +options=addfielddefault(options,'framecolor','black'); +options=addfielddefault(options,'trim','on'); +options=addfielddefault(options,'hardcopy','off'); + +%get fig: +fig=getfieldvalue(options,'figure'); +if ischar(fig), + fig=gcf; +else + figure(fig); + fig=gcf; +end + +%In auto mode, MATLAB prints the figure the same size as it appears on the computer screen, centered on the page +set(fig, 'PaperPositionMode', 'auto'); + +%InvertHardcopy off imposes MATLAB to use the same colors +set(fig, 'InvertHardcopy', getfieldvalue(options,'hardcopy')); + + +%we could have several formats, as a cell array of strings. +formats=format; +if ~iscell(formats), + formats={formats}; +end + +%loop on formats: +for i=1:length(formats), + format=formats{i}; + + %Use higher resolution to anti-alias and use zbuffer to have smooth colors + print(fig, '-zbuffer','-dtiff',['-r' num2str(get(0,'ScreenPixelsPerInch')*getfieldvalue(options,'resolution'))],filename); + + %some trimming involved? + if ~strcmpi(format,'pdf'), + if strcmpi(getfieldvalue(options,'trim'),'on'), + system(['convert -trim ' filename '.tif ' filename '.tif']); + end + end + + %margin? + if ~strcmpi(format,'pdf'), + if strcmpi(getfieldvalue(options,'margin'),'on'), + marginsize=getfieldvalue(options,'marginsize'); + system(['convert -border ' num2str(marginsize) 'x' num2str(marginsize) ' -bordercolor "white" ' filename '.tif ' filename '.tif']); + end + end + + %frame? + if ~strcmpi(format,'pdf'), + if strcmpi(getfieldvalue(options,'frame'),'on'), + framesize=getfieldvalue(options,'framesize'); + framecolor=getfieldvalue(options,'framecolor'); + system(['convert -border ' num2str(framesize) 'x' num2str(framesize) ' -bordercolor "' framecolor '" ' filename '.tif ' filename '.tif']); + end + end + + %convert image to correct format + if ~strcmpi(format,'tiff') & ~strcmpi(format,'tif'), + system(['convert ' filename '.tif ' filename '.' format]); + delete([ filename '.tif']); + end +end Index: sm/trunk-jpl/src/m/model/printmodel.m =================================================================== --- /issm/trunk-jpl/src/m/model/printmodel.m (revision 13005) +++ (revision ) @@ -1,100 +1,0 @@ -function printmodel(filename,format,varargin) -%PRINTMODEL - save an image of current figure -% -% filename: output name of image file (no extension) -% format: image format (ex: 'tiff','jpg','pdf') -% -% List of options to printfmodel: -% -% figure: number of figure to print (default: current figure) -% resolution: use higher resolution to anti-alias (default 2) -% margin: add margin around final image -% marginsize: size of margin around final image (default 5) -% frame: add frame around final image -% framesize: size of frame around final image (default 5) -% framecolor: color of frame around final image (default 'black') -% trim: trim empty space around image (default 'off') -% hardcopy: 'off' to impose MATLAB to use the same colors (default 'off') -% -% Usage: -% printmodel(filename,format,varargin); -% -% Examples: -% printmodel('image','tiff') -% printmodel('image','eps','margin','on','frame','on','hardcopy','on') - - -%get options: -options=pairoptions(varargin{:}); - -%set defaults -options=addfielddefault(options,'figure','gcf'); -options=addfielddefault(options,'format','tiff'); -options=addfielddefault(options,'resolution',2); -options=addfielddefault(options,'margin','on'); -options=addfielddefault(options,'marginsize',5); -options=addfielddefault(options,'frame','on'); -options=addfielddefault(options,'framesize',2); -options=addfielddefault(options,'framecolor','black'); -options=addfielddefault(options,'trim','on'); -options=addfielddefault(options,'hardcopy','off'); - -%get fig: -fig=getfieldvalue(options,'figure'); -if ischar(fig), - fig=gcf; -else - figure(fig); - fig=gcf; -end - -%In auto mode, MATLAB prints the figure the same size as it appears on the computer screen, centered on the page -set(fig, 'PaperPositionMode', 'auto'); - -%InvertHardcopy off imposes MATLAB to use the same colors -set(fig, 'InvertHardcopy', getfieldvalue(options,'hardcopy')); - - -%we could have several formats, as a cell array of strings. -formats=format; -if ~iscell(formats), - formats={formats}; -end - -%loop on formats: -for i=1:length(formats), - format=formats{i}; - - %Use higher resolution to anti-alias and use zbuffer to have smooth colors - print(fig, '-zbuffer','-dtiff',['-r' num2str(get(0,'ScreenPixelsPerInch')*getfieldvalue(options,'resolution'))],filename); - - %some trimming involved? - if ~strcmpi(format,'pdf'), - if strcmpi(getfieldvalue(options,'trim'),'on'), - system(['convert -trim ' filename '.tif ' filename '.tif']); - end - end - - %margin? - if ~strcmpi(format,'pdf'), - if strcmpi(getfieldvalue(options,'margin'),'on'), - marginsize=getfieldvalue(options,'marginsize'); - system(['convert -border ' num2str(marginsize) 'x' num2str(marginsize) ' -bordercolor "white" ' filename '.tif ' filename '.tif']); - end - end - - %frame? - if ~strcmpi(format,'pdf'), - if strcmpi(getfieldvalue(options,'frame'),'on'), - framesize=getfieldvalue(options,'framesize'); - framecolor=getfieldvalue(options,'framecolor'); - system(['convert -border ' num2str(framesize) 'x' num2str(framesize) ' -bordercolor "' framecolor '" ' filename '.tif ' filename '.tif']); - end - end - - %convert image to correct format - if ~strcmpi(format,'tiff') & ~strcmpi(format,'tif'), - system(['convert ' filename '.tif ' filename '.' format]); - delete([ filename '.tif']); - end -end Index: sm/trunk-jpl/src/m/model/qpr.m =================================================================== --- /issm/trunk-jpl/src/m/model/qpr.m (revision 13005) +++ (revision ) @@ -1,7 +1,0 @@ -function qpr(name,format) -%QPR: quick print of a model -% -% Usage: qpr(name,format) - -printmodel(name,format,'trim','on','resolution',1,'margin','on','marginsize',25,'frame','on','framesize',3); - Index: sm/trunk-jpl/src/m/model/setflowequation.m =================================================================== --- /issm/trunk-jpl/src/m/model/setflowequation.m (revision 13005) +++ (revision ) @@ -1,281 +1,0 @@ -function md=setflowequation(md,varargin) -%SETELEMENTSTYPE - associate a solution type to each element -% -% This routine works like plotmodel: it works with an even number of inputs -% 'hutter','macayeal','pattyn','stokes' and 'fill' are the possible options -% that must be followed by the corresponding exp file or flags list -% It can either be a domain file (argus type, .exp extension), or an array of element flags. -% If user wants every element outside the domain to be -% setflowequationd, add '~' to the name of the domain file (ex: '~Pattyn.exp'); -% an empty string '' will be considered as an empty domain -% a string 'all' will be considered as the entire domain -% You can specify the type of coupling, 'penalties' or 'tiling', to use with the input 'coupling' -% -% Usage: -% md=setflowequation(md,varargin) -% -% Example: -% md=setflowequation(md,'pattyn','Pattyn.exp','macayeal',md.mask.elementonfloatingice,'fill','hutter'); -% md=setflowequation(md,'pattyn','Pattyn.exp',fill','hutter','coupling','tiling'); - -%some checks on list of arguments -if ((nargin<2) | (nargout~=1)), - error('setflowequation error message'); -end - -%Find_out what kind of coupling to use -options=pairoptions(varargin{:}); -coupling_method=getfieldvalue(options,'coupling','tiling'); -if (~strcmpi(coupling_method,'tiling') & ~strcmpi(coupling_method,'penalties')), - error('coupling type can only be: tiling or penalties'); -end - -[hutterflag macayealflag pattynflag stokesflag filltype]=recover_areas(md,varargin{:}); - -%Flag the elements that have not been flagged as filltype -if strcmpi(filltype,'hutter'), - hutterflag(find(~(macayealflag | pattynflag)))=1; -elseif strcmpi(filltype,'macayeal'), - macayealflag(find(~(hutterflag | pattynflag | stokesflag)))=1; -elseif strcmpi(filltype,'pattyn'), - pattynflag(find(~(hutterflag | macayealflag | stokesflag)))=1; -end - -%check that each element has at least one flag -if any(hutterflag+ macayealflag+pattynflag+stokesflag==0), - error('setflowequation error message: elements type not assigned, must be specified') -end - -%check that each element has only one flag -if any(hutterflag+ macayealflag+pattynflag+stokesflag>1), - disp('setflowequation warning message: some elements have several types, higher order type is used for them') - hutterflag(find(hutterflag & macayealflag))=0; - hutterflag(find(hutterflag & pattynflag))=0; - macayealflag(find(macayealflag & pattynflag))=0; -end - -%Check that no pattyn or stokes for 2d mesh -if (md.mesh.dimension==2), - if any(stokesflag | pattynflag) - error('setflowequation error message: stokes and pattyn elements not allowed in 2d mesh, extrude it first') - end -end - -%Stokes can only be used alone for now: -if any(stokesflag) &any(hutterflag), - error('setflowequation error message: stokes cannot be used with any other model for now, put stokes everywhere') -end - -%Initialize node fields -nodeonhutter=zeros(md.mesh.numberofvertices,1); -nodeonhutter(md.mesh.elements(find(hutterflag),:))=1; -nodeonmacayeal=zeros(md.mesh.numberofvertices,1); -nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1; -nodeonpattyn=zeros(md.mesh.numberofvertices,1); -nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1; -nodeonstokes=zeros(md.mesh.numberofvertices,1); -noneflag=zeros(md.mesh.numberofelements,1); - -%First modify stokesflag to get rid of elements contrained everywhere (spc + border with pattyn or macayeal) -if any(stokesflag), - fullspcnodes=double((~isnan(md.diagnostic.spcvx)+~isnan(md.diagnostic.spcvy)+~isnan(md.diagnostic.spcvz))==3 | (nodeonpattyn & nodeonstokes)); %find all the nodes on the boundary of the domain without icefront - fullspcelems=double(sum(fullspcnodes(md.mesh.elements),2)==6); %find all the nodes on the boundary of the domain without icefront - stokesflag(find(fullspcelems))=0; - nodeonstokes(md.mesh.elements(find(stokesflag),:))=1; -end - -%Then complete with NoneApproximation or the other model used if there is no stokes -if any(stokesflag), - if any(pattynflag), %fill with pattyn - pattynflag(~stokesflag)=1; - nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1; - elseif any(macayealflag), %fill with macayeal - macayealflag(~stokesflag)=1; - nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1; - else %fill with none - noneflag(find(~stokesflag))=1; - end -end - -%Now take care of the coupling between MacAyeal and Pattyn -md.diagnostic.vertex_pairing=[]; -nodeonmacayealpattyn=zeros(md.mesh.numberofvertices,1); -nodeonpattynstokes=zeros(md.mesh.numberofvertices,1); -nodeonmacayealstokes=zeros(md.mesh.numberofvertices,1); -macayealpattynflag=zeros(md.mesh.numberofelements,1); -macayealstokesflag=zeros(md.mesh.numberofelements,1); -pattynstokesflag=zeros(md.mesh.numberofelements,1); -if strcmpi(coupling_method,'penalties'), - %Create the border nodes between Pattyn and MacAyeal and extrude them - numnodes2d=md.mesh.numberofvertices2d; - numlayers=md.mesh.numberoflayers; - bordernodes2d=find(nodeonpattyn(1:numnodes2d) & nodeonmacayeal(1:numnodes2d)); %Nodes connected to two different types of elements - - %initialize and fill in penalties structure - if ~isnan(bordernodes2d), - penalties=[]; - for i=1:numlayers-1, - penalties=[penalties; [bordernodes2d bordernodes2d+md.mesh.numberofvertices2d*(i)]]; - end - md.diagnostic.vertex_pairing=penalties; - end -elseif strcmpi(coupling_method,'tiling'), - if any(macayealflag) & any(pattynflag), %coupling macayeal pattyn - %Find node at the border - nodeonmacayealpattyn(find(nodeonmacayeal & nodeonpattyn))=1; - %Macayeal elements in contact with this layer become MacAyealPattyn elements - matrixelements=ismember(md.mesh.elements,find(nodeonmacayealpattyn)); - commonelements=sum(matrixelements,2)~=0; - commonelements(find(pattynflag))=0; %only one layer: the elements previously in macayeal - macayealflag(find(commonelements))=0; %these elements are now macayealpattynelements - macayealpattynflag(find(commonelements))=1; - nodeonmacayeal(:)=0; - nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1; - - %rule out elements that don't touch the 2 boundaries - pos=find(macayealpattynflag); - elist=zeros(length(pos),1); - elist = elist + any(sum(nodeonmacayeal(md.mesh.elements(pos,:)),2),2); - elist = elist - any(sum(nodeonpattyn(md.mesh.elements(pos,:)) ,2),2); - pos1=find(elist==1); - macayealflag(pos(pos1))=1; - macayealpattynflag(pos(pos1))=0; - pos2=find(elist==-1); - pattynflag(pos(pos2))=1; - macayealpattynflag(pos(pos2))=0; - - %Recompute nodes associated to these elements - nodeonmacayeal(:)=0; - nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1; - nodeonpattyn(:)=0; - nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1; - nodeonmacayealpattyn(:)=0; - nodeonmacayealpattyn(md.mesh.elements(find(macayealpattynflag),:))=1; - - elseif any(pattynflag) & any(stokesflag), %coupling pattyn stokes - %Find node at the border - nodeonpattynstokes(find(nodeonpattyn & nodeonstokes))=1; - %Stokes elements in contact with this layer become PattynStokes elements - matrixelements=ismember(md.mesh.elements,find(nodeonpattynstokes)); - commonelements=sum(matrixelements,2)~=0; - commonelements(find(pattynflag))=0; %only one layer: the elements previously in macayeal - stokesflag(find(commonelements))=0; %these elements are now macayealpattynelements - pattynstokesflag(find(commonelements))=1; - nodeonstokes=zeros(md.mesh.numberofvertices,1); - nodeonstokes(md.mesh.elements(find(stokesflag),:))=1; - - %rule out elements that don't touch the 2 boundaries - pos=find(pattynstokesflag); - elist=zeros(length(pos),1); - elist = elist + any(sum(nodeonstokes(md.mesh.elements(pos,:)),2),2); - elist = elist - any(sum(nodeonpattyn(md.mesh.elements(pos,:)),2),2); - pos1=find(elist==1); - stokesflag(pos(pos1))=1; - pattynstokesflag(pos(pos1))=0; - pos2=find(elist==-1); - pattynflag(pos(pos2))=1; - pattynstokesflag(pos(pos2))=0; - - %Recompute nodes associated to these elements - nodeonstokes(:)=0; - nodeonstokes(md.mesh.elements(find(stokesflag),:))=1; - nodeonpattyn(:)=0; - nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1; - nodeonpattynstokes(:)=0; - nodeonpattynstokes(md.mesh.elements(find(pattynstokesflag),:))=1; - - elseif any(stokesflag) & any(macayealflag), - %Find node at the border - nodeonmacayealstokes(find(nodeonmacayeal & nodeonstokes))=1; - %Stokes elements in contact with this layer become MacAyealStokes elements - matrixelements=ismember(md.mesh.elements,find(nodeonmacayealstokes)); - commonelements=sum(matrixelements,2)~=0; - commonelements(find(macayealflag))=0; %only one layer: the elements previously in macayeal - stokesflag(find(commonelements))=0; %these elements are now macayealmacayealelements - macayealstokesflag(find(commonelements))=1; - nodeonstokes=zeros(md.mesh.numberofvertices,1); - nodeonstokes(md.mesh.elements(find(stokesflag),:))=1; - - %rule out elements that don't touch the 2 boundaries - pos=find(macayealstokesflag); - elist=zeros(length(pos),1); - elist = elist + any(sum(nodeonmacayeal(md.mesh.elements(pos,:)),2),2); - elist = elist - any(sum(nodeonstokes(md.mesh.elements(pos,:)) ,2),2); - pos1=find(elist==1); - macayealflag(pos(pos1))=1; - macayealstokesflag(pos(pos1))=0; - pos2=find(elist==-1); - stokesflag(pos(pos2))=1; - macayealstokesflag(pos(pos2))=0; - - %Recompute nodes associated to these elements - nodeonmacayeal(:)=0; - nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1; - nodeonstokes(:)=0; - nodeonstokes(md.mesh.elements(find(stokesflag),:))=1; - nodeonmacayealstokes(:)=0; - nodeonmacayealstokes(md.mesh.elements(find(macayealstokesflag),:))=1; - - elseif any(stokesflag) & any(hutterflag), - error('type of coupling not supported yet'); - end -end - -%Create MacaAyealPattynApproximation where needed -md.flowequation.element_equation=zeros(md.mesh.numberofelements,1); -md.flowequation.element_equation(find(noneflag))=0; -md.flowequation.element_equation(find(hutterflag))=1; -md.flowequation.element_equation(find(macayealflag))=2; -md.flowequation.element_equation(find(pattynflag))=3; -md.flowequation.element_equation(find(stokesflag))=4; -md.flowequation.element_equation(find(macayealpattynflag))=5; -md.flowequation.element_equation(find(macayealstokesflag))=6; -md.flowequation.element_equation(find(pattynstokesflag))=7; - -%border -md.flowequation.borderpattyn=nodeonpattyn; -md.flowequation.bordermacayeal=nodeonmacayeal; -md.flowequation.borderstokes=nodeonstokes; - -%Create vertices_type -md.flowequation.vertex_equation=zeros(md.mesh.numberofvertices,1); -pos=find(nodeonhutter); -md.flowequation.vertex_equation(pos)=1; -pos=find(nodeonmacayeal); -md.flowequation.vertex_equation(pos)=2; -pos=find(nodeonpattyn); -md.flowequation.vertex_equation(pos)=3; -pos=find(nodeonhutter); -md.flowequation.vertex_equation(pos)=1; -pos=find(nodeonmacayealpattyn); -md.flowequation.vertex_equation(pos)=5; -pos=find(nodeonstokes); -md.flowequation.vertex_equation(pos)=4; -if any(stokesflag), - pos=find(~nodeonstokes); - if(~any(pattynflag) & ~any(macayealflag)), - md.flowequation.vertex_equation(pos)=0; - end -end -pos=find(nodeonpattynstokes); -md.flowequation.vertex_equation(pos)=7; -pos=find(nodeonmacayealstokes); -md.flowequation.vertex_equation(pos)=6; - -%figure out solution types -md.flowequation.ishutter=double(any(md.flowequation.element_equation==1)); -md.flowequation.ismacayealpattyn=double(any(md.flowequation.element_equation==2 | md.flowequation.element_equation==3)); -md.flowequation.isstokes=double(any(md.flowequation.element_equation==4)); - -return - -%Check that tiling can work: -if any(md.flowequation.bordermacayeal) & any(md.flowequation.borderpattyn) & any(md.flowequation.borderpattyn + md.flowequation.bordermacayeal ~=1), - error('error coupling domain too irregular'); -end -if any(md.flowequation.bordermacayeal) & any(md.flowequation.borderstokes) & any(md.flowequation.borderstokes + md.flowequation.bordermacayeal ~=1), - error('error coupling domain too irregular'); -end -if any(md.flowequation.borderstokes) & any(md.flowequation.borderpattyn) & any(md.flowequation.borderpattyn + md.flowequation.borderstokes~=1), - error('error coupling domain too irregular'); -end Index: sm/trunk-jpl/src/m/model/setflowequation.py =================================================================== --- /issm/trunk-jpl/src/m/model/setflowequation.py (revision 13005) +++ (revision ) @@ -1,278 +1,0 @@ -import numpy -from model import * -from pairoptions import * -from recover_areas import * -from MatlabFuncs import * - -def setflowequation(md,*args): - """ - SETELEMENTSTYPE - associate a solution type to each element - - This routine works like plotmodel: it works with an even number of inputs - 'hutter','macayeal','pattyn','stokes' and 'fill' are the possible options - that must be followed by the corresponding exp file or flags list - It can either be a domain file (argus type, .exp extension), or an array of element flags. - If user wants every element outside the domain to be - setflowequationd, add '~' to the name of the domain file (ex: '~Pattyn.exp'); - an empty string '' will be considered as an empty domain - a string 'all' will be considered as the entire domain - You can specify the type of coupling, 'penalties' or 'tiling', to use with the input 'coupling' - - Usage: - md=setflowequation(md,varargin) - - Example: - md=setflowequation(md,'pattyn','Pattyn.exp','macayeal',md.mask.elementonfloatingice,'fill','hutter'); - md=setflowequation(md,'pattyn','Pattyn.exp',fill','hutter','coupling','tiling'); - """ - - #some checks on list of arguments - if not isinstance(md,model) or not len(args): - raise TypeError("setflowequation error message") - - #Find_out what kind of coupling to use - options=pairoptions(*args) - coupling_method=options.getfieldvalue('coupling','tiling') - if not strcmpi(coupling_method,'tiling') and not strcmpi(coupling_method,'penalties'): - raise TypeError("coupling type can only be: tiling or penalties") - - hutterflag,macayealflag,pattynflag,stokesflag,filltype=recover_areas(md,*args) - - #Flag the elements that have not been flagged as filltype - if strcmpi(filltype,'hutter'): - hutterflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(macayealflag,pattynflag)))]=1 - elif strcmpi(filltype,'macayeal'): - macayealflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(pattynflag,stokesflag))))]=1 - elif strcmpi(filltype,'pattyn'): - pattynflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(macayealflag,stokesflag))))]=1 - - #check that each element has at least one flag - if not numpy.any(hutterflag+macayealflag+pattynflag+stokesflag): - raise TypeError("setflowequation error message: elements type not assigned, must be specified") - - #check that each element has only one flag - if numpy.any(hutterflag+macayealflag+pattynflag+stokesflag>1): - print "setflowequation warning message: some elements have several types, higher order type is used for them" - hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,macayealflag))]=0 - hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,pattynflag))]=0 - macayealflag[numpy.nonzero(numpy.logical_and(macayealflag,pattynflag))]=0 - - #Check that no pattyn or stokes for 2d mesh - if md.mesh.dimension==2: - if numpy.any(numpy.logical_or(stokesflag,pattynflag)): - raise TypeError("setflowequation error message: stokes and pattyn elements not allowed in 2d mesh, extrude it first") - - #Stokes can only be used alone for now: - if numpy.any(stokesflag) and numpy.any(hutterflag): - raise TypeError("setflowequation error message: stokes cannot be used with any other model for now, put stokes everywhere") - - #Initialize node fields - nodeonhutter=numpy.zeros(md.mesh.numberofvertices) - nodeonhutter[md.mesh.elements[numpy.nonzero(hutterflag),:].astype(int)-1]=1 - nodeonmacayeal=numpy.zeros(md.mesh.numberofvertices) - nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1 - nodeonpattyn=numpy.zeros(md.mesh.numberofvertices) - nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1 - nodeonstokes=numpy.zeros(md.mesh.numberofvertices) - noneflag=numpy.zeros(md.mesh.numberofelements) - - #First modify stokesflag to get rid of elements contrained everywhere (spc + border with pattyn or macayeal) - if any(stokesflag): -# fullspcnodes=double((~isnan(md.diagnostic.spcvx)+~isnan(md.diagnostic.spcvy)+~isnan(md.diagnostic.spcvz))==3 | (nodeonpattyn & nodeonstokes)); %find all the nodes on the boundary of the domain without icefront - fullspcnodes=numpy.logical_or(numpy.logical_not(numpy.isnan(md.diagnostic.spcvx))+ \ - numpy.logical_not(numpy.isnan(md.diagnostic.spcvy))+ \ - numpy.logical_not(numpy.isnan(md.diagnostic.spcvz))==3, \ - numpy.logical_and(nodeonpattyn,nodeonstokes)).astype(int) #find all the nodes on the boundary of the domain without icefront -# fullspcelems=double(sum(fullspcnodes(md.mesh.elements),2)==6); %find all the nodes on the boundary of the domain without icefront - fullspcelems=(numpy.sum(fullspcnodes[md.mesh.elements.astype(int)-1],axis=1)==6).astype(int) #find all the nodes on the boundary of the domain without icefront - stokesflag[numpy.nonzero(fullspcelems)]=0 - nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1 - - #Then complete with NoneApproximation or the other model used if there is no stokes - if any(stokesflag): - if any(pattynflag): #fill with pattyn - pattynflag[numpy.logical_not(stokesflag)]=1 - nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1 - elif any(macayealflag): #fill with macayeal - macayealflag[numpy.logical_not(stokesflag)]=1 - nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1 - else: #fill with none - noneflag[numpy.nonzero(numpy.logical_not(stokesflag))]=1 - - #Now take care of the coupling between MacAyeal and Pattyn - md.diagnostic.vertex_pairing=numpy.array([]) - nodeonmacayealpattyn=numpy.zeros(md.mesh.numberofvertices) - nodeonpattynstokes=numpy.zeros(md.mesh.numberofvertices) - nodeonmacayealstokes=numpy.zeros(md.mesh.numberofvertices) - macayealpattynflag=numpy.zeros(md.mesh.numberofelements) - macayealstokesflag=numpy.zeros(md.mesh.numberofelements) - pattynstokesflag=numpy.zeros(md.mesh.numberofelements) - if strcmpi(coupling_method,'penalties'): - #Create the border nodes between Pattyn and MacAyeal and extrude them - numnodes2d=md.mesh.numberofvertices2d - numlayers=md.mesh.numberoflayers - bordernodes2d=numpy.nonzero(numpy.logical_and(nodeonpattyn[1:numnodes2d],nodeonmacayeal[1:numnodes2d])) #Nodes connected to two different types of elements - - #initialize and fill in penalties structure - if numpy.all(numpy.logical_not(numpy.isnan(bordernodes2d))): - penalties=numpy.zeros((0,2)) - for i in xrange(1,numlayers): - penalties=numpy.concatenate((penalties,numpy.concatenate((bordernodes2d,bordernodes2d+md.mesh.numberofvertices2d*(i)),axis=1)),axis=0) - md.diagnostic.vertex_pairing=penalties - - elif strcmpi(coupling_method,'tiling'): - if numpy.any(macayealflag) and numpy.any(pattynflag): #coupling macayeal pattyn - #Find node at the border - nodeonmacayealpattyn[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonpattyn))]=1 - #Macayeal elements in contact with this layer become MacAyealPattyn elements - matrixelements=ismember(md.mesh.elements,numpy.nonzero(nodeonmacayealpattyn)) - commonelements=numpy.sum(matrixelements,axis=1)!=0 - commonelements[numpy.nonzero(pattynflag)]=0 #only one layer: the elements previously in macayeal - macayealflag[numpy.nonzero(commonelements)]=0 #these elements are now macayealpattynelements - macayealpattynflag[numpy.nonzero(commonelements)]=1 - nodeonmacayeal[:]=0 - nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1 - - #rule out elements that don't touch the 2 boundaries - pos=numpy.nonzero(macayealpattynflag) - elist=numpy.zeros(len(pos)) - elist = elist + numpy.any(numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:].astype(int)-1],axis=1),axis=1) - elist = elist - numpy.any(numpy.sum(nodeonpattyn[md.mesh.elements[pos,:].astype(int)-1] ,axis=1),axis=1) - pos1=[i for i,item in enumerate(elist) if item==1] - macayealflag[pos[pos1]]=1 - macayealpattynflag[pos[pos1]]=0 - pos2=[i for i,item in enumerate(elist) if item==-1] - pattynflag[pos[pos2]]=1 - macayealpattynflag[pos[pos2]]=0 - - #Recompute nodes associated to these elements - nodeonmacayeal[:]=0 - nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1 - nodeonpattyn[:]=0 - nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1 - nodeonmacayealpattyn[:]=0 - nodeonmacayealpattyn[md.mesh.elements[numpy.nonzero(macayealpattynflag),:].astype(int)-1]=1 - - elif numpy.any(pattynflag) and numpy.any(stokesflag): #coupling pattyn stokes - #Find node at the border - nodeonpattynstokes[numpy.nonzero(numpy.logical_and(nodeonpattyn,nodeonstokes))]=1 - #Stokes elements in contact with this layer become PattynStokes elements - matrixelements=ismember(md.mesh.elements,numpy.nonzero(nodeonpattynstokes)) - commonelements=numpy.sum(matrixelements,axis=1)!=0 - commonelements[numpy.nonzero(pattynflag)]=0 #only one layer: the elements previously in macayeal - stokesflag[numpy.nonzero(commonelements)]=0 #these elements are now macayealpattynelements - pattynstokesflag[numpy.nonzero(commonelements)]=1 - nodeonstokes=numpy.zeros(md.mesh.numberofvertices) - nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1 - - #rule out elements that don't touch the 2 boundaries - pos=numpy.nonzero(pattynstokesflag) - elist=numpy.zeros(len(pos)) - elist = elist + numpy.any(numpy.sum(nodeonstokes[md.mesh.elements[pos,:].astype(int)-1],axis=1),axis=1) - elist = elist - numpy.any(numpy.sum(nodeonpattyn[md.mesh.elements[pos,:].astype(int)-1],axis=1),axis=1) - pos1=[i for i,item in enumerate(elist) if item==1] - stokesflag[pos[pos1]]=1 - pattynstokesflag[pos[pos1]]=0 - pos2=[i for i,item in enumerate(elist) if item==-1] - pattynflag[pos[pos2]]=1 - pattynstokesflag[pos[pos2]]=0 - - #Recompute nodes associated to these elements - nodeonstokes[:]=0 - nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1 - nodeonpattyn[:]=0 - nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1 - nodeonpattynstokes[:]=0 - nodeonpattynstokes[md.mesh.elements[numpy.nonzero(pattynstokesflag),:].astype(int)-1]=1 - - elif numpy.any(stokesflag) and numpy.any(macayealflag): - #Find node at the border - nodeonmacayealstokes[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonstokes))]=1 - #Stokes elements in contact with this layer become MacAyealStokes elements - matrixelements=ismember(md.mesh.elements,numpy.nonzero(nodeonmacayealstokes)) - commonelements=numpy.sum(matrixelements,axis=1)!=0 - commonelements[numpy.nonzero(macayealflag)]=0 #only one layer: the elements previously in macayeal - stokesflag[numpy.nonzero(commonelements)]=0 #these elements are now macayealmacayealelements - macayealstokesflag[numpy.nonzero(commonelements)]=1 - nodeonstokes=numpy.zeros(md.mesh.numberofvertices) - nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1 - - #rule out elements that don't touch the 2 boundaries - pos=numpy.nonzero(macayealstokesflag) - elist=numpy.zeros(len(pos)) - elist = elist + numpy.any(numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:].astype(int)-1],axis=1),axis=1) - elist = elist - numpy.any(numpy.sum(nodeonstokes[md.mesh.elements[pos,:].astype(int)-1] ,axis=1),axis=1) - pos1=[i for i,item in enumerate(elist) if item==1] - macayealflag[pos[pos1]]=1 - macayealstokesflag[pos[pos1]]=0 - pos2=[i for i,item in enumerate(elist) if item==-1] - stokesflag[pos[pos2]]=1 - macayealstokesflag[pos[pos2]]=0 - - #Recompute nodes associated to these elements - nodeonmacayeal[:]=0 - nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1 - nodeonstokes[:]=0 - nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1 - nodeonmacayealstokes[:]=0 - nodeonmacayealstokes[md.mesh.elements[numpy.nonzero(macayealstokesflag),:].astype(int)-1]=1 - - elif numpy.any(stokesflag) and numpy.any(hutterflag): - raise TypeError("type of coupling not supported yet") - - #Create MacaAyealPattynApproximation where needed - md.flowequation.element_equation=numpy.zeros(md.mesh.numberofelements) - md.flowequation.element_equation[numpy.nonzero(noneflag)]=0 - md.flowequation.element_equation[numpy.nonzero(hutterflag)]=1 - md.flowequation.element_equation[numpy.nonzero(macayealflag)]=2 - md.flowequation.element_equation[numpy.nonzero(pattynflag)]=3 - md.flowequation.element_equation[numpy.nonzero(stokesflag)]=4 - md.flowequation.element_equation[numpy.nonzero(macayealpattynflag)]=5 - md.flowequation.element_equation[numpy.nonzero(macayealstokesflag)]=6 - md.flowequation.element_equation[numpy.nonzero(pattynstokesflag)]=7 - - #border - md.flowequation.borderpattyn=nodeonpattyn - md.flowequation.bordermacayeal=nodeonmacayeal - md.flowequation.borderstokes=nodeonstokes - - #Create vertices_type - md.flowequation.vertex_equation=numpy.zeros(md.mesh.numberofvertices) - pos=numpy.nonzero(nodeonhutter) - md.flowequation.vertex_equation[pos]=1 - pos=numpy.nonzero(nodeonmacayeal) - md.flowequation.vertex_equation[pos]=2 - pos=numpy.nonzero(nodeonpattyn) - md.flowequation.vertex_equation[pos]=3 - pos=numpy.nonzero(nodeonhutter) - md.flowequation.vertex_equation[pos]=1 - pos=numpy.nonzero(nodeonmacayealpattyn) - md.flowequation.vertex_equation[pos]=5 - pos=numpy.nonzero(nodeonstokes) - md.flowequation.vertex_equation[pos]=4 - if numpy.any(stokesflag): - pos=numpy.nonzero(numpy.logical_not(nodeonstokes)) - if not (numpy.any(pattynflag) or numpy.any(macayealflag)): - md.flowequation.vertex_equation[pos]=0 - pos=numpy.nonzero(nodeonpattynstokes) - md.flowequation.vertex_equation[pos]=7 - pos=numpy.nonzero(nodeonmacayealstokes) - md.flowequation.vertex_equation[pos]=6 - - #figure out solution types - md.flowequation.ishutter=float(numpy.any(md.flowequation.element_equation==1)) - md.flowequation.ismacayealpattyn=float(numpy.any(numpy.logical_or(md.flowequation.element_equation==2,md.flowequation.element_equation==3))) - md.flowequation.isstokes=float(numpy.any(md.flowequation.element_equation==4)) - - return md - - #Check that tiling can work: - if numpy.any(md.flowequation.bordermacayeal) and numpy.any(md.flowequation.borderpattyn) and numpy.any(md.flowequation.borderpattyn + md.flowequation.bordermacayeal !=1): - raise TypeError("error coupling domain too irregular") - if numpy.any(md.flowequation.bordermacayeal) and numpy.any(md.flowequation.borderstokes) and numpy.any(md.flowequation.borderstokes + md.flowequation.bordermacayeal !=1): - raise TypeError("error coupling domain too irregular") - if numpy.any(md.flowequation.borderstokes) and numpy.any(md.flowequation.borderpattyn) and numpy.any(md.flowequation.borderpattyn + md.flowequation.borderstokes !=1): - raise TypeError("error coupling domain too irregular") - - return md - Index: sm/trunk-jpl/src/m/model/setmask.m =================================================================== --- /issm/trunk-jpl/src/m/model/setmask.m (revision 13005) +++ (revision ) @@ -1,49 +1,0 @@ -function md=setmask(md,floatingicename,groundedicename) -%SETMASK - establish boundaries between grounded and floating ice. -% -% By default, ice is considered grounded. The contour floatingicename defines nodes -% for which ice is floating. The contour groundedicename defines nodes inside an floatingice, -% that are grounded (ie: ice rises, islands, etc ...) -% All input files are in the Argus format (extension .exp). -% -% Usage: -% md=setmask(md,floatingicename,groundedicename) -% -% Examples: -% md=setmask(md,'all',''); -% md=setmask(md,'Iceshelves.exp','Islands.exp'); - -%some checks on list of arguments -if ((nargin~=3) | (nargout~=1)), - help mask - error('mask error message'); -end - -%Get assigned fields -x=md.mesh.x; -y=md.mesh.y; -elements=md.mesh.elements; - -%Assign elementonfloatingice, elementongroundedice, vertexongroundedice and vertexonfloatingice. Only change at your own peril! This is synchronized heavily with the GroundingLineMigration module. {{{ -elementonfloatingice=FlagElements(md,floatingicename); -elementongroundedice=FlagElements(md,groundedicename); - -%Because groundedice nodes and elements can be included into an floatingice, we need to update. Remember, all the previous -%arrays come from domain outlines that can intersect one another: -elementonfloatingice=double((elementonfloatingice & ~elementongroundedice)); -elementongroundedice=double(~elementonfloatingice); - -%the order here is important. we choose vertexongroundedice as default on the grounding line. -vertexonfloatingice=zeros(md.mesh.numberofvertices,1); -vertexongroundedice=zeros(md.mesh.numberofvertices,1); -vertexongroundedice(md.mesh.elements(find(elementongroundedice),:))=1; -vertexonfloatingice(find(~vertexongroundedice))=1; -%}}} - -%Return: -md.mask.elementonfloatingice=elementonfloatingice; -md.mask.vertexonfloatingice=vertexonfloatingice; -md.mask.elementongroundedice=elementongroundedice; -md.mask.vertexongroundedice=vertexongroundedice; -md.mask.vertexonwater=zeros(md.mesh.numberofvertices,1); -md.mask.elementonwater=zeros(md.mesh.numberofelements,1); Index: sm/trunk-jpl/src/m/model/setmask.py =================================================================== --- /issm/trunk-jpl/src/m/model/setmask.py (revision 13005) +++ (revision ) @@ -1,55 +1,0 @@ -from numpy import * -import FlagElements as fe - -def setmask(md, floatingicename, groundedicename): - #SETMASK - establish boundaries between grounded and floating ice. - # - # By default, ice is considered grounded. The contour floatingicename defines nodes - # for which ice is floating. The contour groundedicename defines nodes inside an floatingice, - # that are grounded (ie: ice rises, islands, etc ...) - # All input files are in the Argus format (extension .exp). - # - # Usage: - # md=setmask(md,floatingicename,groundedicename) - # - # Examples: - # md=setmask(md,'all',''); - # md=setmask(md,'Iceshelves.exp','Islands.exp'); - - #%Get assigned fields - x = md.mesh.x - y = md.mesh.y - elements = md.mesh.elements - - #Assign elementonfloatingice, elementongroundedice, vertexongroundedice and vertexonfloatingice. Only change at your own peril! This is synchronized heavily with the GroundingLineMigration module. {{{ - elementonfloatingice = fe.FlagElements(md, floatingicename) - elementongroundedice = fe.FlagElements(md, groundedicename) - - #Because groundedice nodes and elements can be included into an floatingice, we need to update. Remember, all the previous - #arrays come from domain outlines that can intersect one another: - - elementonfloatingice = logical_and(elementonfloatingice,~elementongroundedice) - elementongroundedice = ~elementonfloatingice - - #the order here is important. we choose vertexongroundedice as default on the grounding line. - vertexonfloatingice = zeros(md.mesh.numberofvertices,'bool') - vertexongroundedice = zeros(md.mesh.numberofvertices,'bool') - - pos=argwhere(elementongroundedice==1) - pos=md.mesh.elements[pos,:]-1 - if pos.size: - vertexongroundedice[pos]=True - - pos=argwhere(~vertexongroundedice) - if pos.size: - vertexonfloatingice[pos]=True; - #%}}} - - #Return: - md.mask.elementonfloatingice = double(elementonfloatingice) - md.mask.vertexonfloatingice = double(vertexonfloatingice) - md.mask.elementongroundedice = double(elementongroundedice) - md.mask.vertexongroundedice = double(vertexongroundedice) - md.mask.vertexonwater = zeros(md.mesh.numberofvertices) - md.mask.elementonwater = zeros(md.mesh.numberofelements) - return md Index: sm/trunk-jpl/src/m/model/setmask2.m =================================================================== --- /issm/trunk-jpl/src/m/model/setmask2.m (revision 13005) +++ (revision ) @@ -1,148 +1,0 @@ -function md=setmask2(md,landname,floatingicename,groundedicename) -%GEOGRAPHY2 - establish land, ice sheet and ice shelf areas in a domains. -% -% Usage: -% md=setmask2(md,landname,floatingicename,groundedicename) -% -% Examples: -% md=setmask2(md,'LandName.exp','Iceshelves.exp','Islands.exp'); - -%Get assigned fields -x=md.mesh.x; -y=md.mesh.y; -elements=md.mesh.elements; - -%recover elements and nodes on land. -if ischar(landname), - [vertexonland,elementonland]=ContourToMesh(elements,x,y,landname,'element and node',2); -elseif isfloat(landname), - if size(landname,1)~=md.mesh.numberofelements, - error('Landname for area must be of same size as number of elements in model'); - end - elementonland=landname; - vertexonland=zeros(md.mesh.numberofvertices,1); - vertexonland(md.mesh.elements(find(elementonland),:))=1; -else - error('Invalid area option option'); -end - -%Now, build the connectivity tables for this mesh. -if size(md.mesh.vertexconnectivity,1)~=md.mesh.numberofvertices, - md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices); -end -if size(md.mesh.elementconnectivity,1)~=md.mesh.numberofelements, - md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity); -end - -%any element with 3 nodes on land should be on land: -elementsonwater=find(~elementonland); -wrongelements=elementsonwater(find(( vertexonland(md.mesh.elements(elementsonwater,1)) + vertexonland(md.mesh.elements(elementsonwater,2)) + vertexonland(md.mesh.elements(elementsonwater,3)) ... - )==3)); -elementonland(wrongelements)=1; - -%any element with its barycentre on land should be on land: (only if landname is an expfile) -if ischar(landname), -weights={[1;1;1],[2;1;1],[1;2;1],[1;1;2]}; - for i=1:length(weights), - xelem=x(md.mesh.elements)*weights{i}/sum(weights{i}); - yelem=y(md.mesh.elements)*weights{i}/sum(weights{i}); - end - baryonland=ContourToNodes(xelem,yelem,landname,1); - pos=find(~baryonland); elementonland(pos)=0; - pos=find(baryonland); elementonland(pos)=1; -end - -%figure out which elements on land are actually in the middle of the ocean! -pos1=find(elementonland); -connectedtoland=md.mesh.elementconnectivity(pos1,:); -pos=find(connectedtoland); connectedtoland(pos)=1-elementonland(connectedtoland(pos)); -connectedtolandsum=sum(connectedtoland,2); -waterelements=pos1(find(connectedtolandsum==3)); -elementonland(waterelements)=0; - -%figure out which elements on water are actually in the middle of the land! -pos1=find(~elementonland); -connectedtowater=md.mesh.elementconnectivity(pos1,:); -pos=find(connectedtowater); connectedtowater(pos)=elementonland(connectedtowater(pos)); -connectedtowatersum=sum(connectedtowater,2); -landelements=pos1(find(connectedtowatersum==3)); -elementonland(landelements)=1; - -%recover arrays of ice shelf nodes and elements, and ice sheet nodes and elements. -elementonfloatingice=FlagElements(md,floatingicename); -elementongroundedice=FlagElements(md,groundedicename); - -%Because groundedice nodes and elements can be included into an floatingice, we need to update. Remember, all the previous -%arrays come from domain outlines that can intersect one another: -vertexonfloatingice=zeros(md.mesh.numberofvertices,1); -vertexongroundedice=zeros(md.mesh.numberofvertices,1); -elementonfloatingice=double((elementonfloatingice & ~elementongroundedice)); -elementongroundedice=double(~elementonfloatingice); -vertexonfloatingice(md.mesh.elements(find(elementonfloatingice),:))=1; -vertexongroundedice(md.mesh.elements(find(elementongroundedice),:))=1; - -%now correct, so that none of the floatingice and groundedice elements and nodes are in the water. -pos=find(~elementonland); -elementonfloatingice(pos)=0; -elementongroundedice(pos)=0; - -pos=find(~vertexonland); -vertexonfloatingice(pos)=0; -vertexongroundedice(pos)=0; - -%create vertexonwater and elementonwater: -vertexonwater=double(~vertexonland); -elementonwater=double(~elementonland); - -%correct for islands: -vertexonfloatingice=double(vertexonfloatingice & ~vertexongroundedice); -elementonfloatingice=double(elementonfloatingice & ~elementongroundedice); - -%now, groundedices are everything except iceshelves and water -vertexongroundedice=double(~vertexonfloatingice & ~vertexonwater); -elementongroundedice=double(~elementonfloatingice & ~elementonwater); - -%Deal with segments on neumann: - -%Get current connectivity -mesh.elementconnectivity=md.mesh.elementconnectivity; - -%put 0 for elements on water -pos=find(mesh.elementconnectivity); -mesh.elementconnectivity(pos)=mesh.elementconnectivity(pos).*(~elementonwater(mesh.elementconnectivity(pos))); - -%put line of ones for elements on water -pos=find(elementonwater); -mesh.elementconnectivity(pos,:)=1;% line of ones for elements on water so they won't be considered - -%resort lines (zeros must be at the last column for findsegments) -mesh.elementconnectivity=sort(mesh.elementconnectivity,2,'descend'); - -%call findsegments to build segment using THIS conectivity -md.mesh.segments=findsegments(md,'mesh.elementconnectivity',mesh.elementconnectivity); - -%some final checks: -%check that no node thinks it's on an ice shelf or ice sheet, and lies actually in the middle of the water. -nodesgrounded=find(~vertexonwater); -lengthconnectivity=size(md.mesh.vertexconnectivity,2); -groundedcounters=md.mesh.vertexconnectivity(nodesgrounded,lengthconnectivity); -groundedconnectivity=md.mesh.vertexconnectivity(nodesgrounded,1:lengthconnectivity-1); -pos=find(groundedconnectivity); -groundedconnectivity(pos)=elementonwater(groundedconnectivity(pos)); -groundedsum=sum(groundedconnectivity,2); -errorflags=find(groundedsum==groundedcounters); -errornodes=nodesgrounded(errorflags); - -vertexonwater(errornodes)=1; -vertexongroundedice(errornodes)=0; -vertexonfloatingice(errornodes)=0; - -%Return: -md.mask.vertexonfloatingice=vertexonfloatingice; -md.mask.elementonfloatingice=elementonfloatingice; -md.mask.vertexonwater=vertexonwater; -md.mask.elementonwater=elementonwater; -md.mask.vertexongroundedice=vertexongroundedice; -md.mask.elementongroundedice=elementongroundedice; - -md.mesh.segmentmarkers(:)=1; Index: /issm/trunk-jpl/src/m/model/solve/parseresultsfromdisk.m =================================================================== --- /issm/trunk-jpl/src/m/model/solve/parseresultsfromdisk.m (revision 13006) +++ /issm/trunk-jpl/src/m/model/solve/parseresultsfromdisk.m (revision 13006) @@ -0,0 +1,212 @@ +function results=parseresultsfromdisk(filename,iosplit) +%PARSERESULTSFROMDISK - ... +% +% Usage: +% results=parseresultsfromdisk(filename,iosplit) + +if iosplit, + results=parseresultsfromdiskiosplit(filename); +else + results=parseresultsfromdiskioserial(filename); +end + +%process patch if necessary +results=MatlabProcessPatch(results); + +function results=parseresultsfromdiskioserial(filename) % {{{ +%PARSERESULTSFROMDISK - ... +% +% Usage: +% results=parseresultsfromdiskioserial(filename) + + +%Open file +fid=fopen(filename,'rb'); +if(fid==-1), + error(['loadresultsfromdisk error message: could not open ',filename,' for binary reading']); +end +results=struct(); + +%Read fields until the end of the file. +result=ReadData(fid); +while ~isempty(result), + %Get time and step + results(result.step).step=result.step; + results(result.step).time=result.time; + + %Add result + if (length(results)>=result.step & isfield(results,result.fieldname) & ~strcmp(result.fieldname,'SolutionType')), + results(result.step).(result.fieldname)=[ results(result.step).(result.fieldname); result.field]; + else + results(result.step).(result.fieldname)=result.field; + end + + %read next result + result=ReadData(fid); + +end + +fclose(fid); +% }}} +function results=parseresultsfromdiskiosplit(filename) % {{{ +%PARSERESULTSFROMDISKIOSPLIT - ... +% +% Usage: +% results=parseresultsfromdiskiosplit(filename) + + +%Open file +fid=fopen(filename,'rb'); +if(fid==-1), + error(['loadresultsfromdisk error message: could not open ',filename,' for binary reading']); +end +results=struct(); + +%if we have done split I/O, ie, we have results that are fragmented across patches, +%do a first pass, and figure out the structure of results +result=ReadDataDimensions(fid); +while ~isempty(result), + + %Get time and step + results(result.step).step=result.step; + results(result.step).time=result.time; + + %Add result + if strcmpi(result.fieldname,'Patch'), + results(result.step).(result.fieldname)=[0 result.N]; + else + results(result.step).(result.fieldname)=NaN; + end + + %read next result + result=ReadDataDimensions(fid); +end + +%do a second pass, and figure out the size of the patches +fseek(fid,0,-1); %rewind +result=ReadDataDimensions(fid); +while ~isempty(result), + + %Add result + if strcmpi(result.fieldname,'Patch'), + patchdimensions=results(result.step).(result.fieldname); + results(result.step).(result.fieldname)=[patchdimensions(1)+result.M result.N]; + end + + %read next result + result=ReadDataDimensions(fid); +end + +%allocate patches +for i=1:length(results), + results(i).Patch=zeros(results(i).Patch(1),results(i).Patch(2)); + results(i).counter=1; %use to index into the patch +end + +%third pass, this time to read the real information +fseek(fid,0,-1); %rewind +result=ReadData(fid); +while ~isempty(result), + + %Get time and step + results(result.step).step=result.step; + results(result.step).time=result.time; + + %Add result + if strcmpi(result.fieldname,'Patch'), + counter=results(result.step).counter; + counter2=counter+size(result.field,1)-1; + results(result.step).(result.fieldname)(counter:counter2,:)=result.field; + + %increment counter: + results(result.step).counter=counter2+1; + else + results(result.step).(result.fieldname)=result.field; + end + + %read next result + result=ReadData(fid); + +end + +%close file +fclose(fid); + % }}} +function result=ReadData(fid) % {{{ +%READDATA - ... +% +% Usage: +% field=ReadData(fid) + +%read field +[length,count]=fread(fid,1,'int'); + +if count==0, + result=struct([]); +else + fieldname=fread(fid,length,'char'); + fieldname=fieldname(1:end-1)'; + fieldname=char(fieldname); + time=fread(fid,1,'double'); + step=fread(fid,1,'int'); + + type=fread(fid,1,'int'); + M=fread(fid,1,'int'); + if type==1, + field=fread(fid,M,'double'); + elseif type==2, + field=fread(fid,M,'char'); + field=char(field(1:end-1)'); + elseif type==3, + N=fread(fid,1,'int'); + field=transpose(fread(fid,[N M],'double')); + else + error(['cannot read data of type ' num2str(type) ]); + end + + result.fieldname=fieldname; + result.time=time; + result.step=step; + result.field=field; +end +% }}} +function result=ReadDataDimensions(fid) % {{{ +%READDATADIMENSIONS - read data dimensions, step and time, but not the data itself. +% +% Usage: +% field=ReadDataDimensions(fid) + + +%read field +[length,count]=fread(fid,1,'int'); + +if count==0, + result=struct([]); +else + fieldname=fread(fid,length,'char'); + fieldname=fieldname(1:end-1)'; + fieldname=char(fieldname); + time=fread(fid,1,'double'); + step=fread(fid,1,'int'); + + type=fread(fid,1,'int'); + M=fread(fid,1,'int'); + N=1; %default + if type==1, + fseek(fid,M*8,0); + elseif type==2, + fseek(fid,M,0); + elseif type==3, + N=fread(fid,1,'int'); + fseek(fid,N*M*8,0); + else + error(['cannot read data of type ' num2str(type) ]); + end + + result.fieldname=fieldname; + result.time=time; + result.step=step; + result.M=M; + result.N=N; +end +% }}} Index: /issm/trunk-jpl/src/m/model/solve/parseresultsfromdisk.py =================================================================== --- /issm/trunk-jpl/src/m/model/solve/parseresultsfromdisk.py (revision 13006) +++ /issm/trunk-jpl/src/m/model/solve/parseresultsfromdisk.py (revision 13006) @@ -0,0 +1,235 @@ +import struct +import numpy +from MatlabFuncs import * +from MatlabProcessPatch import * + +def parseresultsfromdisk(filename,iosplit): + """ + PARSERESULTSFROMDISK - ... + + Usage: + results=parseresultsfromdisk(filename,iosplit) + """ + + if iosplit: + results=parseresultsfromdiskiosplit(filename) + else: + results=parseresultsfromdiskioserial(filename) + + #process patch if necessary + results=MatlabProcessPatch(results) + + return results + +def parseresultsfromdiskioserial(filename): # {{{ + """ + PARSERESULTSFROMDISK - ... + + Usage: + results=parseresultsfromdiskioserial(filename) + """ + + #Open file + try: + fid=open(filename,'rb') + except IOError as e: + raise IOError("loadresultsfromdisk error message: could not open '%s' for binary reading." % filename) + + results={} + + #Read fields until the end of the file. + result=ReadData(fid) + while result: + #Get time and step + if not result['step'] in results: + results[result['step']]={} + results[result['step']]['step']=result['step'] + results[result['step']]['time']=result['time'] + + #Add result + if result['step'] in results and \ + result['fieldname'] in results[result['step']] and \ + not strcmp(result['fieldname'],'SolutionType'): + results[result['step']][result['fieldname']]=numpy.concatenate((results[result['step']][result['fieldname']],result['field']),axis=0) + else: + results[result['step']][result['fieldname']]=result['field'] + + #read next result + result=ReadData(fid) + + fid.close() + + return results + # }}} + +def parseresultsfromdiskiosplit(filename): # {{{ + """ + PARSERESULTSFROMDISKIOSPLIT - ... + + Usage: + results=parseresultsfromdiskiosplit(filename) + """ + + #Open file + try: + fid=open(filename,'rb') + except IOError as e: + raise IOError("loadresultsfromdisk error message: could not open '%s' for binary reading." % filename) + + results={} + + #if we have done split I/O, ie, we have results that are fragmented across patches, + #do a first pass, and figure out the structure of results + result=ReadDataDimensions(fid) + while result: + + #Get time and step + if not result['step'] in results: + results[result['step']]={} + results[result['step']]['step']=result['step'] + results[result['step']]['time']=result['time'] + + #Add result + if strcmpi(result['fieldname'],'Patch'): + results[result['step']][result['fieldname']]=[0,result['N']] + else: + results[result['step']][result['fieldname']]=float('NaN') + + #read next result + result=ReadDataDimensions(fid) + + #do a second pass, and figure out the size of the patches + fid.seek(0) #rewind + result=ReadDataDimensions(fid) + while result: + + #Add result + if strcmpi(result['fieldname'],'Patch'): + patchdimensions=results[result['step']][result['fieldname']] + results[result['step']][result['fieldname']]=[patchdimensions[0]+result['M'],result['N']] + + #read next result + result=ReadDataDimensions(fid) + + #allocate patches + for result in results.itervalues(): + if 'Patch' in result: + result['Patch']=numpy.zeros(shape=(result['Patch'][0],result['Patch'][1]),dtype=float) + result['counter']=0 #use to index into the patch + + #third pass, this time to read the real information + fid.seek(0) #rewind + result=ReadData(fid) + while result: + + #Get time and step + if not result['step'] in results: + results[result['step']]={} + results[result['step']]['step']=result['step'] + results[result['step']]['time']=result['time'] + + #Add result + if strcmpi(result['fieldname'],'Patch'): + counter=results[result['step']]['counter'] + counter2=counter+result['field'].shape[0]-1 + results[result['step']][result['fieldname']][counter:counter2,:]=result['field'] + + #increment counter: + results[result['step']]['counter']=counter2+1 + else: + results[result['step']][result['fieldname']]=result['field'] + + #read next result + result=ReadData(fid) + + #close file + fid.close() + + return results + # }}} + +def ReadData(fid): # {{{ + """ + READDATA - ... + + Usage: + field=ReadData(fid) + """ + + #read field + try: + length=struct.unpack('i',fid.read(struct.calcsize('i')))[0] + + fieldname=struct.unpack('%ds' % length,fid.read(length))[0][:-1] + time=struct.unpack('d',fid.read(struct.calcsize('d')))[0] + step=struct.unpack('i',fid.read(struct.calcsize('i')))[0] + + type=struct.unpack('i',fid.read(struct.calcsize('i')))[0] + M=struct.unpack('i',fid.read(struct.calcsize('i')))[0] + if type==1: + field=numpy.array(struct.unpack('%dd' % M,fid.read(M*struct.calcsize('d'))),dtype=float) + elif type==2: + field=struct.unpack('%ds' % M,fid.read(M))[0][:-1] + elif type==3: + N=struct.unpack('i',fid.read(struct.calcsize('i')))[0] +# field=transpose(fread(fid,[N M],'double')); + field=numpy.zeros(shape=(M,N),dtype=float) + for i in xrange(M): + field[i,:]=struct.unpack('%dd' % N,fid.read(N*struct.calcsize('d'))) + else: + raise TypeError("cannot read data of type %d" % type) + + result={} + result['fieldname']=fieldname + result['time']=time + result['step']=step + result['field']=field + + except struct.error as e: + result={} + + return result + # }}} + +def ReadDataDimensions(fid): # {{{ + """ + READDATADIMENSIONS - read data dimensions, step and time, but not the data itself. + + Usage: + field=ReadDataDimensions(fid) + """ + + #read field + try: + length=struct.unpack('i',fid.read(struct.calcsize('i')))[0] + + fieldname=struct.unpack('%ds' % length,fid.read(length))[0][:-1] + time=struct.unpack('d',fid.read(struct.calcsize('d')))[0] + step=struct.unpack('i',fid.read(struct.calcsize('i')))[0] + + type=struct.unpack('i',fid.read(struct.calcsize('i')))[0] + M=struct.unpack('i',fid.read(struct.calcsize('i')))[0] + N=1 #default + if type==1: + fid.seek(M*8,1) + elif type==2: + fid.seek(M,1) + elif type==3: + N=struct.unpack('i',fid.read(struct.calcsize('i')))[0] + fid.seek(N*M*8,1) + else: + raise TypeError("cannot read data of type %d" % type) + + result={} + result['fieldname']=fieldname + result['time']=time + result['step']=step + result['M']=M + result['N']=N + + except struct.error as e: + result={} + + return result + # }}} + Index: /issm/trunk-jpl/src/m/model/solve/waitonlock.m =================================================================== --- /issm/trunk-jpl/src/m/model/solve/waitonlock.m (revision 13006) +++ /issm/trunk-jpl/src/m/model/solve/waitonlock.m (revision 13006) @@ -0,0 +1,58 @@ +function flag=waitonlock(md,executionpath,login,port) +%WAITONLOCK - wait for a file +% +% This routine will return when a file named 'filename' is written to disk. +% If the time limit given in input is exceeded, return 0 +% +% Usage: +% flag=waitonlock(md,executionpath) + +%Get filename (lock file) and options +executionpath=md.cluster.executionpath; +cluster=md.cluster.name; +login=md.cluster.login; +port=md.cluster.port; +timelimit=md.settings.waitonlock; +filename=[executionpath '/' md.private.runtimename '/' md.miscellaneous.name '.lock']; + +%waitonlock will work if the lock is on the same machine only: +if ~strcmpi(oshostname(),cluster), + + disp('solution launched on remote cluster. log in to detect job completion.'); + choice=input('Is the job successfully completed? (y/n)','s'); + if ~strcmp(choice,'y'), + disp('Results not loaded... exiting'); + flag=0; + else + flag=1; + end + +%job is running on the same machine +else + + if ismember('interactive',properties(md.cluster)) & md.cluster.interactive + %We are in interactive mode, no need to check for job completion + flag=1; + return; + end + %initialize time and file presence test flag + time=0; ispresent=0; + disp(['waiting for ' filename ' hold on... (Ctrl+C to exit)']) + + %loop till file .lock exist or time is up + while (ispresent==0 & timetimelimit), + disp('Time limit exceeded. Increase md.settings.waitonlock'); + disp('The results must be loaded manually with md=loadresultsfromcluster(md).'); + error(['waitonlock error message: time limit exceeded']); + flag=0; + else + flag=1; + end +end Index: /issm/trunk-jpl/src/m/model/solve/waitonlock.py =================================================================== --- /issm/trunk-jpl/src/m/model/solve/waitonlock.py (revision 13006) +++ /issm/trunk-jpl/src/m/model/solve/waitonlock.py (revision 13006) @@ -0,0 +1,64 @@ +import os +import socket +import time +from MatlabFuncs import * + +def waitonlock(md,executionpath,login,port): + """ + WAITONLOCK - wait for a file + + This routine will return when a file named 'filename' is written to disk. + If the time limit given in input is exceeded, return 0 + + Usage: + flag=waitonlock(md,executionpath) + """ + + #Get filename (lock file) and options + executionpath=md.cluster.executionpath + cluster=md.cluster.name + login=md.cluster.login + port=md.cluster.port + timelimit=md.settings.waitonlock + filename=os.path.join(executionpath,md.private.runtimename,md.miscellaneous.name+'.lock') + + #waitonlock will work if the lock is on the same machine only: + if not strcmpi(socket.gethostname().lower().split('.')[0],cluster): + + print 'solution launched on remote cluster. log in to detect job completion.' + choice=raw_input('Is the job successfully completed? (y/n) ') + if not strcmp(choice,'y'): + print 'Results not loaded... exiting' + flag=0 + else: + flag=1 + + #job is running on the same machine + else: + + if 'interactive' in vars(md.cluster) and md.cluster.interactive: + #We are in interactive mode, no need to check for job completion + flag=1 + return flag + #initialize time and file presence test flag + etime=0 + ispresent=0 + print "waiting for '%s' hold on... (Ctrl+C to exit)" % filename + + #loop till file .lock exist or time is up + while ispresent==0 and etimetimelimit: + print 'Time limit exceeded. Increase md.settings.waitonlock' + print 'The results must be loaded manually with md=loadresultsfromcluster(md).' + raise RuntimeError('waitonlock error message: time limit exceeded.') + flag=0 + else: + flag=1 + + return flag + Index: sm/trunk-jpl/src/m/model/waitonlock.m =================================================================== --- /issm/trunk-jpl/src/m/model/waitonlock.m (revision 13005) +++ (revision ) @@ -1,58 +1,0 @@ -function flag=waitonlock(md,executionpath,login,port) -%WAITONLOCK - wait for a file -% -% This routine will return when a file named 'filename' is written to disk. -% If the time limit given in input is exceeded, return 0 -% -% Usage: -% flag=waitonlock(md,executionpath) - -%Get filename (lock file) and options -executionpath=md.cluster.executionpath; -cluster=md.cluster.name; -login=md.cluster.login; -port=md.cluster.port; -timelimit=md.settings.waitonlock; -filename=[executionpath '/' md.private.runtimename '/' md.miscellaneous.name '.lock']; - -%waitonlock will work if the lock is on the same machine only: -if ~strcmpi(oshostname(),cluster), - - disp('solution launched on remote cluster. log in to detect job completion.'); - choice=input('Is the job successfully completed? (y/n)','s'); - if ~strcmp(choice,'y'), - disp('Results not loaded... exiting'); - flag=0; - else - flag=1; - end - -%job is running on the same machine -else - - if ismember('interactive',properties(md.cluster)) & md.cluster.interactive - %We are in interactive mode, no need to check for job completion - flag=1; - return; - end - %initialize time and file presence test flag - time=0; ispresent=0; - disp(['waiting for ' filename ' hold on... (Ctrl+C to exit)']) - - %loop till file .lock exist or time is up - while (ispresent==0 & timetimelimit), - disp('Time limit exceeded. Increase md.settings.waitonlock'); - disp('The results must be loaded manually with md=loadresultsfromcluster(md).'); - error(['waitonlock error message: time limit exceeded']); - flag=0; - else - flag=1; - end -end Index: sm/trunk-jpl/src/m/model/waitonlock.py =================================================================== --- /issm/trunk-jpl/src/m/model/waitonlock.py (revision 13005) +++ (revision ) @@ -1,64 +1,0 @@ -import os -import socket -import time -from MatlabFuncs import * - -def waitonlock(md,executionpath,login,port): - """ - WAITONLOCK - wait for a file - - This routine will return when a file named 'filename' is written to disk. - If the time limit given in input is exceeded, return 0 - - Usage: - flag=waitonlock(md,executionpath) - """ - - #Get filename (lock file) and options - executionpath=md.cluster.executionpath - cluster=md.cluster.name - login=md.cluster.login - port=md.cluster.port - timelimit=md.settings.waitonlock - filename=os.path.join(executionpath,md.private.runtimename,md.miscellaneous.name+'.lock') - - #waitonlock will work if the lock is on the same machine only: - if not strcmpi(socket.gethostname().lower().split('.')[0],cluster): - - print 'solution launched on remote cluster. log in to detect job completion.' - choice=raw_input('Is the job successfully completed? (y/n) ') - if not strcmp(choice,'y'): - print 'Results not loaded... exiting' - flag=0 - else: - flag=1 - - #job is running on the same machine - else: - - if 'interactive' in vars(md.cluster) and md.cluster.interactive: - #We are in interactive mode, no need to check for job completion - flag=1 - return flag - #initialize time and file presence test flag - etime=0 - ispresent=0 - print "waiting for '%s' hold on... (Ctrl+C to exit)" % filename - - #loop till file .lock exist or time is up - while ispresent==0 and etimetimelimit: - print 'Time limit exceeded. Increase md.settings.waitonlock' - print 'The results must be loaded manually with md=loadresultsfromcluster(md).' - raise RuntimeError('waitonlock error message: time limit exceeded.') - flag=0 - else: - flag=1 - - return flag -