source: issm/trunk-jpl/src/m/model/regionaltransient2d.m@ 12618

function md2=regionaltransient2d(md1,area,hmin,hmax,err,stepres)
%regionaltransient2d - extract a model according to an Argus contour or flag list and remesh
%               at new resolution res
%
%   This routine extracts a submodel from a bigger model with respect to a given contour
%   md must be followed by the corresponding exp domain file (argus type, .exp extension). 
%   The model will be remeshed at high rsolution hmin and low resolution hmax.  The ice 
%   boundary velocities will be spc'd to the transient velocities at saved transient steps
%   at the resolution optionally provided for stepres.  A stepres of 2 means that you wish
%   to skip every other saved transient step.  This is useful when extracting a long transient.
%
%   Usage:
%      md2=regionaltransient2d(md1,area,hmin,hmax,err);
%
%   Examples:
%      md2=regionaltransient2d(md,'Domain.exp',500,10000,[15 250]);
%      md2=regionaltransient2d(md,'Domain.exp',3000,15000,[10 300],2);
%
%   See also: MODELEXTRACT, EXTRUDE, COLLAPSE

%some checks
if ((nargin~=5) & (nargin~=6)),
        help regionaltransient2d 
        error('regionaltransient2d error message: bad usage');
end

%get check option
if (nargin==5),
        stepres=1;
end

%take every fields from model
mde=modelextract(md1,area);
mde.private.bamg=[];
mde.mesh.extractedvertices=nan;
mde.mesh.extractedelements=nan;

%remesh
md2=bamg(mde,'hmin',hmin,'hmax',hmax,'field',[mde.inversion.vel_obs mde.geometry.surface],'splitcorner',1,'KeepVertices',0,'err',err);
md2=setmask(md2,'','');

%automatically modify fields

        %loop over model fields
        model_fields=fields(md1);
        for i=1:length(model_fields),

                %get field
                field=md1.(model_fields{i});
                fieldsize=size(field);

                %copy field, interpolated to new mesh
                if isobject(field), %recursive call
                        object_fields=fields(md1.(model_fields{i}));
                        fname=['(model_fields{i}).(object_fields{j})'];
                else
                        object_fields=field;
                        fname=['(model_fields{i})'];
                end
                for j=1:length(object_fields),
                        %get field
                        field=eval(['md2.' fname]);
                        fieldsize=size(field);

                        %size = number of nodes * n
                        for n=1:fieldsize(2)
                                if fieldsize(1)==mde.mesh.numberofvertices
                                        if(sum(field(:,n) ~= field(1,n)) == 0)
                                                eval(['md2.' fname '(1:md2.mesh.numberofvertices,n)=field(1,n)*ones(md2.mesh.numberofvertices,1);']);
                                        else
                                                eval(['md2.' fname '(1:md2.mesh.numberofvertices,n)=InterpFromMeshToMesh2d(mde.mesh.elements,mde.mesh.x,mde.mesh.y,field(:,n),md2.mesh.x,md2.mesh.y);']);
                                        end
                                        eval(['md2.' fname '(:,n)=md2.' fname '(1:md2.mesh.numberofvertices,n);']);
                                elseif fieldsize(1)==mde.mesh.numberofvertices+1
                                        if(sum(field(1:end-1,n) ~= field(1,n)) == 0)
                                                eval(['md2.' fname '(1:md2.mesh.numberofvertices+1,n)=[field(1,n)*ones(md2.mesh.numberofvertices,1); field(end,n)];']);
                                        else
                                                eval(['md2.' fname '(1:md2.mesh.numberofvertices+1,n)=[InterpFromMeshToMesh2d(mde.mesh.elements,mde.mesh.x,mde.mesh.y,field(1:end-1,n),md2.mesh.x,md2.mesh.y); field(end,n)];']);
                                        end
                                        eval(['md2.' fname '(:,n)=md2.' fname '(1:md2.mesh.numberofvertices+1,n)']);
                                        %size = number of elements * n
                                elseif fieldsize(1)==mde.mesh.numberofelements
                                        if(sum(field(1:end-1,n) ~= field(1,n)) == 0)
                                                eval(['md2.' fname '(1:md2.mesh.numberofelements,n)=field(1,n)*ones(md2.mesh.numberofelements,1);']);
                                        else
                                                eval(['md2.' fname '(1:md2.mesh.numberofelements,n)=InterpFromMeshToMesh2d(mde.mesh.elements,mde.mesh.x,mde.mesh.y,field(:,n),md2.mesh.x,md2.mesh.y);']);
                                        end
                                        eval(['md2.' fname '(:,n)=md2.' fname '(1:md2.mesh.numberofelements,n);']);
                                end
                        end
                end
        end

        %Read transient velocities and thickness, looping through only the populated times
        spcx=[];
        spcy=[];
        spct=[];
        steps=[];
        nsteps=length(md1.results.TransientSolution);
        count=0;
        numElements=arrayfun(@(x) numel(x.step), md1.results.TransientSolution);
        for t=find(numElements==1)
                if ~isempty(md1.results.TransientSolution(t).Vel) & mod(count,stepres)==0,
                        vx=PatchToVec(md1.results.TransientSolution(t).Vx);
                        vy=PatchToVec(md1.results.TransientSolution(t).Vy);
                        thickness=PatchToVec(md1.results.TransientSolution(t).Thickness);
                        spcx=[spcx InterpFromMeshToMesh2d(md1.mesh.elements,md1.mesh.x,md1.mesh.y,vx,md2.mesh.x,md2.mesh.y)];
                        spcy=[spcy InterpFromMeshToMesh2d(md1.mesh.elements,md1.mesh.x,md1.mesh.y,vy,md2.mesh.x,md2.mesh.y)];
                        spct=[spct InterpFromMeshToMesh2d(md1.mesh.elements,md1.mesh.x,md1.mesh.y,thickness,md2.mesh.x,md2.mesh.y)];
                        steps=[steps t*md1.timestepping.time_step];
                end
                count=count+1;
        end

        %As long as there are recorded time steps, spc the boundaries with velocities
        if nsteps > 0
                md2.diagnostic.spcvx=md2.diagnostic.spcvx*ones(1,size(spcx,2));
                md2.diagnostic.spcvy=md2.diagnostic.spcvy*ones(1,size(spcy,2));
                md2.diagnostic.spcvz=md2.diagnostic.spcvz*ones(1,size(spcx,2));
                md2.prognostic.spcthickness=md2.prognostic.spcthickness*ones(1,size(spct,2));
                md2.diagnostic.spcvx(find(md2.mesh.vertexonboundary),:)=spcx(find(md2.mesh.vertexonboundary),:);
                md2.diagnostic.spcvy(find(md2.mesh.vertexonboundary),:)=spcy(find(md2.mesh.vertexonboundary),:);
                md2.diagnostic.spcvz(find(md2.mesh.vertexonboundary),:)=0;
                md2.prognostic.spcthickness(find(md2.mesh.vertexonboundary),:)=spct(find(md2.mesh.vertexonboundary),:);
                md2.diagnostic.spcvx=[md2.diagnostic.spcvx; steps];
                md2.diagnostic.spcvy=[md2.diagnostic.spcvy; steps];
                md2.diagnostic.spcvz=[md2.diagnostic.spcvz; steps];
                md2.prognostic.spcthickness=[md2.prognostic.spcthickness; steps];
        end

        %Diagnostic.  Don't spc the icefront vertices.
        if ~isnan(md2.diagnostic.icefront)
                md1s=modelextract(md1,area);
                %md2.diagnostic.icefront=[md2.mesh.segments 2];
                e2=md2.mesh.segments(:,end);
                e1=md1s.mesh.segments(:,end);

                pload = nan*ones(size(md1s.mesh.elements,1),1);
                pload(md1s.diagnostic.icefront(:,end-1))=md1s.diagnostic.icefront(:,end);

                x2=mean(md2.mesh.x(md2.mesh.elements(e2,:)),2);
      y2=mean(md2.mesh.y(md2.mesh.elements(e2,:)),2);
                x1=mean(md1s.mesh.x(md1s.mesh.elements),2);
      y1=mean(md1s.mesh.y(md1s.mesh.elements),2);

                pload2=griddata(x1,y1,pload,x2,y2,'nearest');
                md2.diagnostic.icefront=[md2.mesh.segments(~isnan(pload2),:) pload2(~isnan(pload2))];
                md2.diagnostic.spcvx(unique(md2.diagnostic.icefront(:,1:2)),:)=nan;
                md2.diagnostic.spcvy(unique(md2.diagnostic.icefront(:,1:2)),:)=nan;
                md2.diagnostic.spcvz(unique(md2.diagnostic.icefront(:,1:2)),:)=nan;
                md2.prognostic.spcthickness(unique(md2.diagnostic.icefront(:,1:2)),:)=nan;
        end

        %Clear results fields
        if isstruct(md1.results),
                md2.results=[];
        end