source: issm/trunk/src/m/model/modelextract.m@ 11237

function md2=modelextract(md1,area)
%modelextract - extract a model according to an Argus contour or flag list
%
%   This routine extracts a submodel from a bigger model with respect to a given contour
%   md 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 
%   modelextract2d, 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
%   add an argument 0 if you do not want the elements to be checked (faster)
%
%   Usage:
%      md2=modelextract(md1,area);
%
%   Examples:
%      md2=modelextract(md,'Domain.exp');
%      md2=modelextract(md,md.mask.elementonfloatingice);
%
%   See also: EXTRUDE, COLLAPSE

%some checks
if ((nargin~=2) | (nargout~=1)),
        help modelextract
        error('modelextract error message: bad usage');
end

%get check option
if (nargin==3 & varargin{1}==0),
        checkoutline=0;
else
        checkoutline=1;
end

%get elements that are inside area
flag_elem=FlagElements(md1,area);
if ~any(flag_elem),
        error('extracted model is empty');
end

%kick out all elements with 3 dirichlets
spc_elem=find(~flag_elem);
spc_node=sort(unique(md1.mesh.elements(spc_elem,:)));
flag=ones(md1.mesh.numberofvertices,1);
flag(spc_node)=0;
pos=find(sum(flag(md1.mesh.elements),2)==0);
flag_elem(pos)=0;

%extracted elements and nodes lists
pos_elem=find(flag_elem);
pos_node=sort(unique(md1.mesh.elements(pos_elem,:)));

%keep track of some fields
numberofvertices1=md1.mesh.numberofvertices;
numberofelements1=md1.mesh.numberofelements;
numberofvertices2=length(pos_node);
numberofelements2=length(pos_elem);
flag_node=zeros(numberofvertices1,1);
flag_node(pos_node)=1;

%Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements)
Pelem=zeros(numberofelements1,1);
Pelem(pos_elem)=[1:numberofelements2]';
Pnode=zeros(numberofvertices1,1);
Pnode(pos_node)=[1:numberofvertices2]';

%renumber the elements (some node won't exist anymore)
elements_1=md1.mesh.elements;
elements_2=elements_1(pos_elem,:);
elements_2(:,1)=Pnode(elements_2(:,1));
elements_2(:,2)=Pnode(elements_2(:,2));
elements_2(:,3)=Pnode(elements_2(:,3));
if md1.mesh.dimension==3,
        elements_2(:,4)=Pnode(elements_2(:,4));
        elements_2(:,5)=Pnode(elements_2(:,5));
        elements_2(:,6)=Pnode(elements_2(:,6));
end

%OK, now create the new model !

        %take every fields from model
        md2=md1;

        %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);
                if isobject(field), %recursive call
                        object_fields=fields(md1.(model_fields{i}));
                        for j=1:length(object_fields),
                                %get field
                                field=md1.(model_fields{i}).(object_fields{j});
                                fieldsize=size(field);
                                %size = number of nodes * n
                                if fieldsize(1)==numberofvertices1
                                        md2.(model_fields{i}).(object_fields{j})=field(pos_node,:);
                                elseif (fieldsize(1)==numberofvertices1+1)
                                        md2.(model_fields{i}).(object_fields{j})=[field(pos_node,:); field(end,:)];
                                        %size = number of elements * n
                                elseif fieldsize(1)==numberofelements1
                                        md2.(model_fields{i}).(object_fields{j})=field(pos_elem,:);
                                end
                        end
                else
                        %size = number of nodes * n
                        if fieldsize(1)==numberofvertices1
                                md2.(model_fields{i})=field(pos_node,:);
                        elseif (fieldsize(1)==numberofvertices1+1)
                                md2.(model_fields{i})=[field(pos_node,:); field(end,:)];
                                %size = number of elements * n
                        elseif fieldsize(1)==numberofelements1
                                md2.(model_fields{i})=field(pos_elem,:);
                        end
                end
        end

        %modify some specific fields

        %Mesh
        md2.mesh.numberofelements=numberofelements2;
        md2.mesh.numberofvertices=numberofvertices2;
        md2.mesh.elements=elements_2;

        %mesh.uppervertex mesh.lowervertex
        if md1.mesh.dimension==3
                md2.mesh.uppervertex=md1.mesh.uppervertex(pos_node);
                pos=find(~isnan(md2.mesh.uppervertex));
                md2.mesh.uppervertex(pos)=Pnode(md2.mesh.uppervertex(pos));

                md2.mesh.lowervertex=md1.mesh.lowervertex(pos_node);
                pos=find(~isnan(md2.mesh.lowervertex));
                md2.mesh.lowervertex(pos)=Pnode(md2.mesh.lowervertex(pos));

                md2.mesh.upperelements=md1.mesh.upperelements(pos_elem);
                pos=find(~isnan(md2.mesh.upperelements));
                md2.mesh.upperelements(pos)=Pelem(md2.mesh.upperelements(pos));

                md2.mesh.lowerelements=md1.mesh.lowerelements(pos_elem);
                pos=find(~isnan(md2.mesh.lowerelements));
                md2.mesh.lowerelements(pos)=Pelem(md2.mesh.lowerelements(pos));
        end

        %Initial 2d mesh 
        if md1.mesh.dimension==3
                flag_elem_2d=flag_elem(1:md1.mesh.numberofelements2d);
                pos_elem_2d=find(flag_elem_2d);
                flag_node_2d=flag_node(1:md1.mesh.numberofvertices2d);
                pos_node_2d=find(flag_node_2d);

                md2.mesh.numberofelements2d=length(pos_elem_2d);
                md2.mesh.numberofvertices2d=length(pos_node_2d);
                md2.mesh.elements2d=md1.mesh.elements2d(pos_elem_2d,:);
                md2.mesh.elements2d(:,1)=Pnode(md2.mesh.elements2d(:,1));
                md2.mesh.elements2d(:,2)=Pnode(md2.mesh.elements2d(:,2));
                md2.mesh.elements2d(:,3)=Pnode(md2.mesh.elements2d(:,3));

                md2.mesh.x2d=md1.mesh.x(pos_node_2d);
                md2.mesh.y2d=md1.mesh.y(pos_node_2d);
        end

        %Edges
        if size(md2.mesh.edges,2)>1, %do not use ~isnan because there are some NaNs...
                %renumber first two columns
                pos=find(~isnan(md2.mesh.edges(:,4)));
                md2.mesh.edges(:  ,1)=Pnode(md2.mesh.edges(:,1)); 
                md2.mesh.edges(:  ,2)=Pnode(md2.mesh.edges(:,2)); 
                md2.mesh.edges(:  ,3)=Pelem(md2.mesh.edges(:,3));
                md2.mesh.edges(pos,4)=Pelem(md2.mesh.edges(pos,4));
                %remove edges when the 2 vertices are not in the domain.
                md2.mesh.edges=md2.mesh.edges(find(md2.mesh.edges(:,1) & md2.mesh.edges(:,2)),:);
                %Replace all zeros by NaN in the last two columns;
                pos=find(md2.mesh.edges(:,3)==0);
                md2.mesh.edges(pos,3)=NaN;
                pos=find(md2.mesh.edges(:,4)==0);
                md2.mesh.edges(pos,4)=NaN;
                %Invert NaN of the third column with last column (Also invert first two columns!!)
                pos=find(isnan(md2.mesh.edges(:,3)));
                md2.mesh.edges(pos,3)=md2.mesh.edges(pos,4);
                md2.mesh.edges(pos,4)=NaN;
                values=md2.mesh.edges(pos,2);
                md2.mesh.edges(pos,2)=md2.mesh.edges(pos,1);
                md2.mesh.edges(pos,1)=values;
                %Finally remove edges that do not belong to any element
                pos=find(isnan(md2.mesh.edges(:,3)) & isnan(md2.mesh.edges(:,4)));
                md2.mesh.edges(pos,:)=[];
        end

        %Penalties
        if ~isnan(md2.diagnostic.vertex_pairing),
                for i=1:size(md1.diagnostic.vertex_pairing,1);
                        md2.diagnostic.vertex_pairing(i,:)=Pnode(md1.diagnostic.vertex_pairing(i,:));
                end
                md2.diagnostic.vertex_pairing=md2.diagnostic.vertex_pairing(find(md2.diagnostic.vertex_pairing(:,1)),:);
        end
        if ~isnan(md2.prognostic.vertex_pairing),
                for i=1:size(md1.prognostic.vertex_pairing,1);
                        md2.prognostic.vertex_pairing(i,:)=Pnode(md1.prognostic.vertex_pairing(i,:));
                end
                md2.prognostic.vertex_pairing=md2.prognostic.vertex_pairing(find(md2.prognostic.vertex_pairing(:,1)),:);
        end

        %recreate segments
        if md1.mesh.dimension==2
                md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices);
                md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity);
                md2.mesh.segments=contourenvelope(md2);
                md2.mesh.vertexonboundary=zeros(numberofvertices2,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1;
        end

        %Boundary conditions: Dirichlets on new boundary
        %Catch the elements that have not been extracted
        orphans_elem=find(~flag_elem);
        orphans_node=unique(md1.mesh.elements(orphans_elem,:))';
        %Figure out which node are on the boundary between md2 and md1
        nodestoflag1=intersect(orphans_node,pos_node);
        nodestoflag2=Pnode(nodestoflag1);
        if numel(md1.diagnostic.spcvx)>1 & numel(md1.diagnostic.spcvy)>2 & numel(md1.diagnostic.spcvz)>2,
                if numel(md1.inversion.vx_obs)>1 & numel(md1.inversion.vy_obs)>1
                        md2.diagnostic.spcvx(nodestoflag2)=md2.inversion.vx_obs(nodestoflag2); 
                        md2.diagnostic.spcvy(nodestoflag2)=md2.inversion.vy_obs(nodestoflag2);
                else
                        md2.diagnostic.spcvx(nodestoflag2)=NaN;
                        md2.diagnostic.spcvy(nodestoflag2)=NaN;
                        disp(' ')
                        disp('!! modelextract warning: spc values should be checked !!')
                        disp(' ')
                end
                %put 0 for vz
                md2.diagnostic.spcvz(nodestoflag2)=0;
        end
        if ~isnan(md1.thermal.spctemperature),
                md2.thermal.spctemperature(nodestoflag2,1)=1;
        end

        %Diagnostic
        if ~isnan(md2.diagnostic.icefront)
                md2.diagnostic.icefront(:,1)=Pnode(md1.diagnostic.icefront(:,1)); 
                md2.diagnostic.icefront(:,2)=Pnode(md1.diagnostic.icefront(:,2)); 
                md2.diagnostic.icefront(:,end-1)=Pelem(md1.diagnostic.icefront(:,end-1));
                if md1.mesh.dimension==3
                        md2.diagnostic.icefront(:,3)=Pnode(md1.diagnostic.icefront(:,3)); 
                        md2.diagnostic.icefront(:,4)=Pnode(md1.diagnostic.icefront(:,4)); 
                end
                md2.diagnostic.icefront=md2.diagnostic.icefront(find(md2.diagnostic.icefront(:,1) & md2.diagnostic.icefront(:,2) & md2.diagnostic.icefront(:,end)),:);
        end

        %Results fields
        if isstruct(md1.results),
                md2.results=struct();
                solutionfields=fields(md1.results);
                for i=1:length(solutionfields),
                        %get subfields
                        solutionsubfields=fields(md1.results.(solutionfields{i}));
                        for j=1:length(solutionsubfields),
                                field=md1.results.(solutionfields{i}).(solutionsubfields{j});
                                if length(field)==numberofvertices1,
                                        md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_node);
                                elseif length(field)==numberofelements1,
                                        md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_elem);
                                else
                                        md2.results.(solutionfields{i}).(solutionsubfields{j})=field;
                                end
                        end
                end
        end

%Keep track of pos_node and pos_elem
md2.mesh.extractedvertices=pos_node;
md2.mesh.extractedelements=pos_elem;