source: issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/modelmerge.m@ 21368

function md=modelmerge(md1,md2)
%MODELMERGE  - merge two models by merging their meshes
%
%   Usage:
%      md=modelmerge(md1,md2);

        if nargin~=2 
                help meshconvert
                error('meshconvert error message: bad usage');
        end
        
        md=model();

        %first ,copy md1 mesh into md.mesh to initialize: 
        md.mesh=md1.mesh;

        %some initializatoin: 
        elements1=md1.mesh.elements;
        x1=md1.mesh.x;
        y1=md1.mesh.y;
        nods1=md1.mesh.numberofvertices;
        nel1=md1.mesh.numberofelements;

        elements2=md2.mesh.elements;
        x2=md2.mesh.x;
        y2=md2.mesh.y;
        nods2=md2.mesh.numberofvertices;
        nel2=md2.mesh.numberofelements;
        segs2=md2.mesh.segments;

        %offset elements2 by nods1: 
        elements2=elements2+nods1;

        tolerance=1e-5;
        %go into the vertices on boundary of mesh 1, and figure out which ones are common to mesh2: 
        verticesonboundary=find(md1.mesh.vertexonboundary); 
        for i=1:length(verticesonboundary),
                node1=verticesonboundary(i); xnode1=x1(node1); ynode1=y1(node1);
                %is there another node with these coordinates in mesh2? 
                ind=find(sqrt(((x2-xnode1).^2+(y2-ynode1).^2))<tolerance);
                if ~isempty(ind),
                        x2(ind)=NaN;
                        y2(ind)=NaN;
                        pos=find(elements2==(ind+nods1)); elements2(pos)=node1;
                end
        end

        %go through elements2 and drop counter on each vertex that is above the x2 and y2 vertices being dropped: 
        while( ~isempty(find(isnan(x2)))),
                for i=1:length(x2),
                        if isnan(x2(i)),
                                pos=find(elements2>(i+nods1));
                                elements2(pos)=elements2(pos)-1;
                                x2(i)=[];
                                y2(i)=[];
                                break;
                        end
                end
        end

        %merge elements: 
        elements=[elements1;elements2];

        %merge vertices: 
        x=[x1;x2]; 
        y=[y1;y2];

        %output: 
        md.mesh.x=x;
        md.mesh.y=y;
        md.mesh.elements=elements;
        md.mesh.numberofvertices=length(x);
        md.mesh.numberofelements=size(elements,1);

        %connectivities: 
        md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
        md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);

        %find segments: 
        md.mesh.segments=findsegments(md);

        %vertex on boundary: 
        md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1);
        md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;

        %some checks: 
        if max(md.mesh.elements)>md.mesh.numberofvertices, 
                error('issue in modelmerge, one of the element ids is > number of vertices!');
        end