Changeset 6381


Ignore:
Timestamp:
10/21/10 16:00:46 (15 years ago)
Author:
jschierm
Message:

Working version of kml_mesh_write.m, but needs a bunch of cleaning up.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/m/kml/kml_mesh_write.m

    r6313 r6381  
    262262    [xseg,yseg]=flagedges(md.elements,md.x,md.y,md.part);
    263263    fprintf(fid,'    <Folder>\n');
    264     fprintf(fid,'      <name>Segments</name>\n');
     264    fprintf(fid,'      <name>Partition Segments</name>\n');
    265265    fprintf(fid,'      <visibility>1</visibility>\n');
    266266    fprintf(fid,'      <description>Partitions=%d, Segments=%d</description>\n',...
    … …  
    356356        irow=unique(irow);
    357357        elem=md.elements(irow,:);
    358         nodecon=NodeConnectivity(elem,md.numberofgrids);
     358        nodecon=nodeconnectivity(elem,md.numberofgrids);
    359359        [edgeper,elemper,iloop]=edgeperimeter(elem,nodecon);
    360360        iloop(end+1)=size(edgeper,1)+1;
    … …  
    415415end
    416416
     417%  write folder for partition edges
     418
     419if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
     420    md.npart
     421    fprintf(fid,'    <Folder>\n');
     422    fprintf(fid,'      <name>Partition Edges</name>\n');
     423    fprintf(fid,'      <visibility>1</visibility>\n');
     424    fprintf(fid,'      <description>Partitions=%d, Nodes=%d</description>\n',...
     425        md.npart,md.numberofgrids);
     426
     427%  write each partition loop as linestrings
     428
     429    disp(['Writing ' num2str(md.npart) ' partitions as KML linestrings.']);
     430    epart=md.part(md.elements)+1;
     431    if exist('ndata','var') || exist('edata','var')
     432        pdata=zeros(1,md.npart);
     433        pdata(:)=NaN;
     434    end
     435
     436%  loop over each partition
     437
     438    for k=1:md.npart
     439        disp(['partition k=' int2str(k)])
     440       
     441%  for each partition, find all the included elements and determine the
     442%  perimeter (including those shared by another partition)
     443
     444        [icol,irow]=find(epart'==k);
     445        irow=unique(irow);
     446        elemp=md.elements(irow,:);
     447        epartp=epart(irow,:);
     448        nodeconp=nodeconnectivity(elemp,md.numberofgrids);
     449        [edgeadjp]=edgeadjacency(elemp,nodeconp);
     450        [edgeper,elemper,iloop]=edgeperimeter(elemp,nodeconp,edgeadjp);
     451        iloop(end+1)=size(edgeper,1)+1;
     452       
     453%  determine the data to be used for the colors (if any)
     454
     455        if exist('ndata','var')
     456            pdata(k)=ndata(find(md.part+1==k,1));
     457        elseif exist('edata','var')
     458            for i=1:size(epartp,1)
     459                if isempty(find(epart(i,:)~=k,1))
     460                    pdata(k)=edata(irow(i));
     461                    break
     462                end
     463            end
     464            if isnan(pdata(k))
     465                warning('Data for Partition %d is not defined.\n',k)
     466            end
     467        end
     468       
     469%  loop over each loop of the perimeter for the given partition
     470
     471        for i=1:length(iloop)-1
     472            fprintf(fid,'      <Placemark>\n');
     473            if (length(iloop)-1 > 1)
     474                fprintf(fid,'        <name>Partition %d, Loop %d</name>\n',k,i);
     475            else
     476                fprintf(fid,'        <name>Partition %d</name>\n',k);
     477            end
     478            fprintf(fid,'        <visibility>1</visibility>\n');
     479            if exist('pdata','var')
     480                fprintf(fid,'        <description>Partition data: %g</description>\n',pdata(k));
     481                imap = fix((pdata(k)-cmin)/(cmax-cmin)*size(cmap,1))+1;
     482                if     (imap >= 1) && (imap <= size(cmap,1))
     483                    fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',imap);
     484                elseif (pdata(k) == cmax)
     485                    fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',size(cmap,1));
     486                else
     487                    fprintf(fid,'        <styleUrl>#BlackLineEmptyPoly</styleUrl>\n');
     488                end
     489            else
     490                fprintf(fid,'        <styleUrl>#BlackLineRandomPoly</styleUrl>\n');
     491            end
     492%            fprintf(fid,'        <Polygon>\n');
     493%            fprintf(fid,'          <extrude>1</extrude>\n');
     494%            fprintf(fid,'          <altitudeMode>relativeToGround</altitudeMode>\n');
     495%            fprintf(fid,'          <outerBoundaryIs>\n');
     496%            fprintf(fid,'            <LinearRing>\n');
     497        fprintf(fid,'        <LineString>\n');
     498        fprintf(fid,'          <extrude>1</extrude>\n');
     499        fprintf(fid,'          <tessellate>1</tessellate>\n');
     500        fprintf(fid,'          <altitudeMode>relativeToGround</altitudeMode>\n');
     501            fprintf(fid,'              <coordinates>\n');
     502            elast=0;
     503            nlast=0;
     504            slast=0;
     505            lat=[];
     506            long=[];
     507           
     508%  loop over the element edges on the loop of the partition
     509
     510            j=iloop(i);
     511            while (j < iloop(i+1))
     512%  find which edge of element is referenced in perimeter list
     513                for l=1:size(elemp,2)
     514                    if ((elemp(elemper(j),l)          == edgeper(j,1)) && ...
     515                        (elemp(elemper(j),mod(l,3)+1) == edgeper(j,2))) || ...
     516                       ((elemp(elemper(j),l)          == edgeper(j,2)) && ...
     517                        (elemp(elemper(j),mod(l,3)+1) == edgeper(j,1)))
     518                        jedge=l;
     519                        break
     520                    end
     521                end
     522
     523%  check if element edge connects nodes in partition
     524                if (epartp(elemper(j),jedge)          == k) && ...
     525                   (epartp(elemper(j),mod(jedge,3)+1) == k)
     526%  write out specified element edge
     527                    disp(['segment j=' int2str(j) ' unshared edge ' int2str(edgeper(j,1)) ' to ' int2str(edgeper(j,2)) ' on side ' int2str(jedge) ' from element ' int2str(elemper(j)) ' written.'])
     528                    if ~elast
     529                        [lat(end+1),long(end+1)]=mapxy(md.x(edgeper(j,1)),md.y(edgeper(j,1)),'s');
     530                        fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
     531                    end
     532                    [lat(end+1),long(end+1)]=mapxy(md.x(edgeper(j,2)),md.y(edgeper(j,2)),'s');
     533                    fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
     534                    elast=elemper(j);
     535                    nlast=edgeper(j,2);
     536                    slast=0;
     537                    j=j+1;
     538                   
     539%  element not entirely within partition, so figure out boundary
     540                else
     541                    disp(['segment j=' int2str(j) ' from element ' int2str(elemper(j)) ' shared by other partitions.'])
     542                    ielem=elemper(j);
     543
     544%  follow partition boundary through elements not wholly in partition
     545%  (may include elements not in perimeter list)
     546
     547                    while 1
     548%  if first edge, figure out direction from perimeter direction
     549                        if ~nlast && ~slast
     550                            nlast=find(elemp(ielem,:)==edgeper(j,1));
     551                            nnext=find(elemp(ielem,:)==edgeper(j,2));
     552                            if     (nlast+nnext == 3)
     553                                slast=1;
     554                            elseif (nlast+nnext == 5)
     555                                slast=2;
     556                            elseif (nlast+nnext == 4)
     557                                slast=3;
     558                            end
     559                            if     (nnext+(6-nlast-nnext) == 3)
     560                                snext=1;
     561                            elseif (nnext+(6-nlast-nnext) == 5)
     562                                snext=2;
     563                            elseif (nnext+(6-nlast-nnext) == 4)
     564                                snext=3;
     565                            end
     566
     567%  find how many nodes of current element are in current partition
     568%  (1 or 2, not 3)
     569                            ipart=find(epartp(ielem,:)==k);
     570%  two nodes are in current partition, so cut off other node
     571                            if (length(ipart) == 2)
     572                                switch 6-sum(ipart)
     573                                    case nlast
     574                                        slast=6-slast-snext;
     575                                        nlast=0;
     576                                    case nnext
     577                                        if (epartp(ielem,nnext) == k)
     578                                            nlast=nnext;
     579                                        end
     580                                    otherwise
     581                                        slast=6-slast-snext;
     582                                        nlast=0;
     583                                end
     584%  one node is in current partition
     585                            else
     586%  all different, so cut through centroid
     587                                if (epartp(ielem,1) ~= epartp(ielem,2)) && ...
     588                                   (epartp(ielem,2) ~= epartp(ielem,3)) && ...
     589                                   (epartp(ielem,3) ~= epartp(ielem,1))
     590                                    switch ipart
     591                                        case {nlast,nnext}
     592                                            if (epartp(ielem,nnext) == k)
     593                                                nlast=nnext;
     594                                            end
     595                                        otherwise
     596                                            slast=6-slast-snext;
     597                                            nlast=0;
     598                                    end
     599                                else
     600%  other two are in the same partition, so cut them off
     601                                    switch ipart
     602                                        case nlast
     603                                            if (epartp(ielem,nnext) == k)
     604                                                nlast=nnext;
     605                                            end
     606                                        case nnext
     607                                            slast=snext;
     608                                            nlast=0;
     609                                        otherwise
     610                                            slast=6-slast-snext;
     611                                            nlast=0;
     612                                    end
     613                                end
     614                            end
     615
     616%  last edge exited last element at node
     617                            if nlast
     618%  write out first node of first edge for half-edge to midpoint
     619                                disp(['segment j=' int2str(j) ' unshared half edge from node ' int2str(elemp(ielem,nlast)) ' (node ' int2str(nlast) ') on side ' int2str(slast) ' from element ' int2str(ielem) ' written.'])
     620                                [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,nlast))),...
     621                                                               (md.y(elemp(ielem,nlast))),'s');
     622                                fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
     623                            end
     624                            nlast=0;
     625                           
     626                            if ~elast
     627%  write out midpoint of first edge
     628                                [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,slast))...
     629                                                               +md.x(elemp(ielem,mod(slast,3)+1)))/2.,...
     630                                                               (md.y(elemp(ielem,slast))...
     631                                                               +md.y(elemp(ielem,mod(slast,3)+1)))/2.,'s');
     632                                fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
     633                            end
     634                        end
     635
     636%  last edge exited last element at node
     637                        if nlast
     638                            if elast
     639%  find where last node on previous element occurs on current element
     640                                nlast=find(elemp(ielem,:)==nlast,1);
     641                            end
     642%  half-edge occurs on unshared edge from current node (unique unless mesh
     643%  is only attached at node)
     644                            switch nlast
     645                                case 1
     646                                    if ~edgeadjp(ielem,1)
     647                                        nnext=2;
     648                                        slast=1;
     649                                    else
     650                                        nnext=3;
     651                                        slast=3;
     652                                    end
     653                                case 2
     654                                    if ~edgeadjp(ielem,2)
     655                                        nnext=3;
     656                                        slast=2;
     657                                    else
     658                                        nnext=1;
     659                                        slast=1;
     660                                    end
     661                                case 3
     662                                    if ~edgeadjp(ielem,3)
     663                                        nnext=1;
     664                                        slast=3;
     665                                    else
     666                                        nnext=2;
     667                                        slast=2;
     668                                    end
     669                            end
     670%  write out half-edge from current node to midpoint of unshared edge
     671                            disp(['segment j=' int2str(j) ' unshared half edge from node ' int2str(elemp(ielem,nlast)) ' (node ' int2str(nlast) ') on side ' int2str(slast) ' from element ' int2str(ielem) ' written.'])
     672                            [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,nlast))...
     673                                                           +md.x(elemp(ielem,nnext)))/2.,...
     674                                                           (md.y(elemp(ielem,nlast))...
     675                                                           +md.y(elemp(ielem,nnext)))/2.,'s');
     676                            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
     677                            nlast=0;
     678
     679%  last edge exited last element at midpoint of side
     680                        elseif slast
     681                            if elast
     682%  find where last edge on previous element occurs on current element
     683                                slast=find(edgeadjp(ielem,:)==elast,1);
     684                            end
     685                        end
     686
     687%  find how many nodes of current element are in current partition
     688%  (1 or 2, not 3)
     689                        ipart=find(epartp(ielem,:)==k);
     690                        if (length(ipart) == 2)
     691%  two nodes are in current partition, so cut off other node
     692                            switch 6-sum(ipart)
     693                                case 1
     694                                    snext=3+1-slast;
     695                                case 2
     696                                    snext=1+2-slast;
     697                                case 3
     698                                    snext=2+3-slast;
     699                            end
     700                        else
     701                            if (epartp(ielem,1) ~= epartp(ielem,2)) && ...
     702                               (epartp(ielem,2) ~= epartp(ielem,3)) && ...
     703                               (epartp(ielem,3) ~= epartp(ielem,1))
     704%  all different, so write out centroid
     705                                disp(['element ielem=' int2str(ielem) ' centroid written.'])
     706                                [lat(end+1),long(end+1)]=mapxy(sum(md.x(elemp(ielem,:)))/3.,...
     707                                                               sum(md.y(elemp(ielem,:)))/3.,'s');
     708                                fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
     709                            end
     710%  one node is in current partition, so cut off other two nodes
     711                            switch ipart
     712                                case 1
     713                                    snext=3+1-slast;
     714                                case 2
     715                                    snext=1+2-slast;
     716                                case 3
     717                                    snext=2+3-slast;
     718                            end
     719                        end
     720%  write out midpoint of opposite edge
     721                        disp(['segment j=' int2str(j) ' internal edge from side ' int2str(slast) ' to side ' int2str(snext) ' from element ' int2str(ielem) ' written.'])
     722                        [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,snext))...
     723                                                       +md.x(elemp(ielem,mod(snext,3)+1)))/2.,...
     724                                                       (md.y(elemp(ielem,snext))...
     725                                                       +md.y(elemp(ielem,mod(snext,3)+1)))/2.,'s');
     726                        fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
     727                        elast=ielem;
     728                        nlast=0;
     729                        slast=snext;
     730%  find adjacent element to opposite edge
     731                        ielem=edgeadjp(elast,slast);
     732%  if opposite edge is unshared, find it in edge perimeter list
     733                        if ~ielem
     734                            jlast=find(elemper(j:end)==elast)+j-1;
     735                            j=0;
     736                            for l=1:length(jlast)
     737                                if ((elemp(elast,slast)          == edgeper(jlast(l),1)) && ...
     738                                    (elemp(elast,mod(slast,3)+1) == edgeper(jlast(l),2))) || ...
     739                                   ((elemp(elast,slast)          == edgeper(jlast(l),2)) && ...
     740                                    (elemp(elast,mod(slast,3)+1) == edgeper(jlast(l),1)))
     741                                    j=jlast(l);
     742                                    break
     743                                end
     744                            end
     745                            if ~j
     746                                j=iloop(i+1)-1;
     747                            end
     748%  write out half-edge from midpoint of unshared edge to node
     749                            if (epartp(elast,slast) == k)
     750                                nnext=slast;
     751                            else
     752                                nnext=mod(slast,3)+1;
     753                            end
     754                            disp(['segment j=' int2str(j) ' unshared half edge on side ' int2str(slast) ' to node ' int2str(elemp(elast,nnext)) ' (node ' int2str(nnext) ') from element ' int2str(elast) ' written.'])
     755                            [lat(end+1),long(end+1)]=mapxy(md.x(elemp(elast,nnext)),...
     756                                                           md.y(elemp(elast,nnext)),'s');
     757                            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
     758                            break
     759%  if not unshared, advance perimeter list and watch for end
     760                        else
     761                            if (elast == elemper(j))
     762                                if (j+1 < iloop(i+1)) && ...
     763                                   ~isempty(find(elemper(j+1:end)~=elast,1))
     764                                    j=j+find(elemper(j+1:end)~=elast,1);
     765                                else
     766                                    break
     767                                end
     768                            end
     769                        end
     770                    end
     771                    j=j+1;
     772                end
     773            end
     774%            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(iloop(i)),lat(iloop(i)),alt);
     775
     776            fprintf(fid,'              </coordinates>\n');
     777        fprintf(fid,'        </LineString>\n');
     778%            fprintf(fid,'            </LinearRing>\n');
     779%            fprintf(fid,'          </outerBoundaryIs>\n');
     780%            fprintf(fid,'        </Polygon>\n');
     781            fprintf(fid,'      </Placemark>\n');
     782        end
     783    end
     784    fprintf(fid,'    </Folder>\n');
     785end
     786
    417787%  write trailer data
    418788
Note: See TracChangeset for help on using the changeset viewer.