source: issm/trunk-jpl/src/m/mesh/meshband.m@ 22938

function band=meshband(mh,outerdomain,varargin)

        %process options: 
        options=pairoptions(varargin{:});

        %some checks on the mesh: 
        if (isempty(mh.x) | isempty(mh.y)  ) 
                error('meshband error message: mesh has one of the following empty: ''x'',''y''');
        end

        %give ourselves a unique temporary directory: 
        temproot=tempname; mkdir(temproot);

        %align external segments on internal mesh: 
        mh.segments=alignsegments(mh.segments);

        %figure out domain outline: 
        meshtodomain(mh,[temproot '/innerdomain.exp']);

        %create domain outine from inner and outer domain
        inner=expread([temproot '/innerdomain.exp']);
        if inner.closed==0,
                inner.nods=inner.nods+1;
                inner.x(end+1)=inner.x(1);
                inner.y(end+1)=inner.y(1);
                inner.closed=1;
                if getfieldvalue(options,'invert',0),
                        inner.x=flipud(inner.x);
                        inner.y=flipud(inner.y);
                end
        end

        [path,name,ext]=fileparts(outerdomain);
        if strcmpi(ext,'.shp'),
                outer=shpread(outerdomain);
        else
                outer=expread(outerdomain);
        end
        if outer.closed==0,
                outer.nods=outer.nods+1;
                outer.x(end+1)=outer.x(1);
                outer.y(end+1)=outer.y(1);
                outer.closed=1;
        end

        domain(1).x=outer.x; 
        domain(1).y=outer.y; 
        domain(1).density=outer.density; 
        domain(1).nods=outer.nods; 
        domain(1).closed=outer.closed;
        
        
        domain(2).x=inner.x; 
        domain(2).y=inner.y; 
        domain(2).density=inner.density; 
        domain(2).nods=inner.nods; 
        domain(2).closed=inner.closed;


        expwrite(domain,[temproot '/Band.exp']);
        
        if getfieldvalue(options,'plot',0),
                figure(1),clf,hold on,axis image;
                expdisp([temproot '/Band.exp'],'linestyle','k-*');
        end

        %mesh: 
        md=bamg(model(),'domain',[temproot '/Band.exp'],'MaxCornerAngle',1e-15,'gradation',10000); band=md.mesh; clear md;

        if getfieldvalue(options,'plot',0),
                figure(2),clf,trisurf(band.elements,band.x,band.y,band.x),view(2),shading faceted;
                hold on,expdisp([temproot '/Band.exp'],'linestyle','k-*');
        end

        %check that the domain vertices = the number of segments: 
        if abs(length(band.segments) - (length(domain(1).x)+length(domain(2).x)-2))>=2,
                disp(sprintf('band mesh not consistent: %i!=%i+%i\n',length(band.segments), length(domain(1).x), length(domain(2).x)));

                figure(3),clf,expdisp([temproot '/Band.exp'],'linestyle','r*');
                hold on; 
                for i=1:length(band.segments),
                        i1=band.segments(i,1); i2=band.segments(i,2);
                        plot([band.x(i1) band.x(i2)],[band.y(i1) band.y(i2)],'k*-');
                        plot([band.x(i1)+ band.x(i2)]/2,[band.y(i1)+ band.y(i2)]/2,'g*');
                end
                %error('band mesh not consistent');
        end

        %erase temporary directory: 
        system(['rm -rf ' temproot]);

end