tricontour.m

Go to the documentation of this file.

/* (Autoinserted by mtoc++)
 * This source code has been filtered by the mtoc++ executable,
 * which generates code that can be processed by the doxygen documentation tool.
 *
 * On the other hand, it can neither be interpreted by MATLAB, nor can it be compiled with a C++ compiler.
 * Except for the comments, the function bodies of your M-file functions are untouched.
 * Consequently, the FILTER_SOURCE_FILES doxygen switch (default in our Doxyfile.template) will produce
 * attached source files that are highly readable by humans.
 *
 * Additionally, links in the doxygen generated documentation to the source code of functions and class members refer to
 * the correct locations in the source code browser.
 * However, the line numbers most likely do not correspond to the line numbers in the original MATLAB source files.
 */

function [cout , hout ] = tricontour(ax,double t,p,Hn,N,varargin) {

/*  Contouring for functions defined on triangular meshes
 *
 *   TRICONTOUR(p,t,F,N)
 *
 * Draws contours of the surface F, where F is defined on the triangulation
 * [p,t]. These inputs define the xy co-ordinates of the points and their
 * connectivity:
 *
 *   P = [x1,y1; x2,y2; etc],            - xy co-ordinates of nodes in the 
 *                                         triangulation
 *   T = [n11,n12,n13; n21,n23,n23; etc] - node numbers in each triangle
 *
 * The last input N defines the contouring levels. There are several
 * options:
 *
 *   N scalar - N number of equally spaced contours will be drawn
 *   N vector - Draws contours at the levels specified in N
 *
 * A special call with a two element N where both elements are equal draws a
 * single contour at that level.
 *
 *   [C,H] = TRICONTOUR(...)
 *
 * This syntax can be used to pass the contour matrix C and the contour
 * handels H to clabel by adding clabel(c,h) or clabel(c) after the call to
 * TRICONTOUR.
 *
 * TRICONTOUR can also return 3D contours similar to CONTOUR3 by adding
 * view(3) after the call to TRICONTOUR.
 *
 * Type "contourdemo" for some examples.
 *
 * See also, CONTOUR, CLABEL */

/*  This function does NOT interpolate back onto a Cartesian grid, but
 * instead uses the triangulation directly.
 *
 * If your going to use this inside a loop with the same [p,t] a good
 * modification is to make the connectivity "mkcon" once outside the loop
 * because usually about 50% of the time is spent in "mkcon".
 *
 * Darren Engwirda - 2005 (d_engwirda@hotmail.com)
 * Updated 15/05/2006 */


/*  I/O checking */
if nargin<5
    error(" Incorrect number of inputs ")
end
if nargout>2
    error(" Incorrect number of outputs ")
end

/*  Error checking */
if (size(p,2)~=2) || (size(t,2)~=3) || (size(Hn,2)~=1)
    error(" Incorrect input dimensions ")
end
if size(p,1)~=size(Hn,1)
    error(" F and p must be the same length ")
end
if (max(t(:))>size(p,1)) || (min(t(:))<=0)
    error(" t is not a valid triangulation of p ")
end
if (size(N,1)>1) && (size(N,2)>1)
    error(" N cannot be a matrix ")
end

/*  Make mesh connectivity data structures (edge based pointers) */
[e,eINt,e2t] = mkcon(p,t);

numt = size(t,1);       /*  Num triangles */

nume = size(e,1);       /*  Num edges */




/* ==========================================================================
 *                Quadratic interpolation to centroids
 *========================================================================== */

/*  Nodes */
t1 = t(:,1); t2 = t(:,2); t3 = t(:,3);

/*  FORM FEM GRADIENTS
 * Evaluate centroidal gradients (piecewise-linear interpolants) */
x23 = p(t2,1)-p(t3,1);  y23 = p(t2,2)-p(t3,2);
x21 = p(t2,1)-p(t1,1);  y21 = p(t2,2)-p(t1,2);

/*  Centroidal values */
Htx = (y23.*Hn(t1) + (y21-y23).*Hn(t2) - y21.*Hn(t3)) ./ (x23.*y21-x21.*y23);
Hty = (x23.*Hn(t1) + (x21-x23).*Hn(t2) - x21.*Hn(t3)) ./ (y23.*x21-y21.*x23);

/*  Form nodal gradients.
 * Take the average of the neighbouring centroidal values */
Hnx = 0*Hn; Hny = Hnx; count = Hnx;
for k = 1:numt
    /*  Nodes */
    n1 = t1(k); n2 = t2(k); n3 = t3(k);
    /*  Current values */
    Hx = Htx(k); Hy = Hty(k);
    /*  Average to n1 */
    Hnx(n1)   = Hnx(n1)+Hx;
    Hny(n1)   = Hny(n1)+Hy;
    count(n1) = count(n1)+1;
    /*  Average to n2 */
    Hnx(n2)   = Hnx(n2)+Hx;
    Hny(n2)   = Hny(n2)+Hy;
    count(n2) = count(n2)+1;
    /*  Average to n3 */
    Hnx(n3)   = Hnx(n3)+Hx;
    Hny(n3)   = Hny(n3)+Hy;
    count(n3) = count(n3)+1;
end
Hnx = Hnx./count;
Hny = Hny./count;

/*  Centroids [x,y] */
pt = (p(t1,:)+p(t2,:)+p(t3,:))/3;

/*  Take unweighted average of the linear extrapolation from nodes to centroids */
Ht = ( Hn(t1) + (pt(:,1)-p(t1,1)).*Hnx(t1) + (pt(:,2)-p(t1,2)).*Hny(t1) + ...
       Hn(t2) + (pt(:,1)-p(t2,1)).*Hnx(t2) + (pt(:,2)-p(t2,2)).*Hny(t2) + ...
       Hn(t3) + (pt(:,1)-p(t3,1)).*Hnx(t3) + (pt(:,2)-p(t3,2)).*Hny(t3) )/3;



/*  DEAL WITH CONTOURING LEVELS */
if length(N)==1
    lev = linspace(max(Ht),min(Ht),N+1);
    num = N;
else
    if (length(N)==2) && (N(1)==N(2))
        lev = N(1);
        num = 1;
    else
        lev = sort(N);
        num = length(N);
        lev = lev(num:-1:1);
    end
end

/*  MAIN LOOP */
c   = [];
h   = [];
in  = false(numt,1);
vec = 1:numt;
old = in;
for v = 1:num       /*  Loop over contouring levels */


    /*  Find centroid values >= current level */
    i     = vec(Ht>=lev(v));
    i     = i(~old(i));         /*  Don't need to check triangles from higher levels */

    in(i) = true;

    /*  Locate boundary edges in group */
    bnd  = [i; i; i];       /*  Just to alloc */

    next = 1;
    for k = 1:length(i)
        ct    = i(k);
        count = 0;
        for q = 1:3     /*  Loop through edges in ct */

            ce = eINt(ct,q);
            if ~in(e2t(ce,1)) || ((e2t(ce,2)>0)&&~in(e2t(ce,2)))    
                bnd(next) = ce;     /*  Found bnd edge */

                next      = next+1;
            else
                count = count+1;    /*  Count number of non-bnd edges in ct */

            end
        end
        if count==3                 /*  If 3 non-bnd edges ct must be in middle of group */

            old(ct) = true;         /*  & doesn't need to be checked for the next level */

        end
    end
    numb = next-1; bnd(next:end) = [];

    /*  Skip to next lev if empty */
    if numb==0
        continue
    end

    /*  Place nodes approximately on contours by interpolating across bnd
     * edges     */
    t1  = e2t(bnd,1);
    t2  = e2t(bnd,2);
    ok  = t2>0;

    /*  Get two points for interpolation. Always use t1 centroid and 
     * use t2 centroid for internal edges and bnd midpoint for boundary 
     * edges */

    /*  1st point is always t1 centroid */
    H1 = Ht(t1);                                                /*  Centroid value */

    p1 = ( p(t(t1,1),:)+p(t(t1,2),:)+p(t(t1,3),:) )/3;          /*  Centroid [x,y] */


    /*  2nd point is either t2 centroid or bnd edge midpoint */
    i1        = t2(ok);                                         /*  Temp indexing */

    i2        = bnd(~ok);
    H2        = H1;
    H2(ok)    = Ht(i1);                                         /*  Centroid values internally */

    H2(~ok)   = ( Hn(e(i2,1))+Hn(e(i2,2)) )/2;                  /*  Edge values at boundary */

    p2        = p1;
    p2(ok,:)  = ( p(t(i1,1),:)+p(t(i1,2),:)+p(t(i1,3),:) )/3;   /*  Centroid [x,y] internally */

    p2(~ok,:) = ( p(e(i2,1),:)+p(e(i2,2),:) )/2;                /*  Edge [x,y] at boundary */


    /*  Linear interpolation */
    r     = (lev(v)-H1)./(H2-H1);
    penew = p1 + [r,r].*(p2-p1);

    /*  Do a temp connection between adjusted node & endpoint nodes in
     * ce so that the connectivity between neighbouring adjusted nodes
     * can be determined */
    vecb    = (1:numb)^t;
    m       = 2*vecb-1;
    c1      = 0*m;
    c2      = 0*m;
    c1(m)   = e(bnd,1);
    c1(m+1) = e(bnd,2);
    c2(m)   = vecb;
    c2(m+1) = vecb;

    /*  Sort connectivity to place connected edges in sucessive rows */
    [c1,i] = sort(c1); c2 = c2(i);

    /*  Connect adjacent adjusted nodes */
    k    = 1;
    next = 1;
    while k<(2*numb)
        if c1(k)==c1(k+1)
            c1(next) = c2(k);
            c2(next) = c2(k+1);
            next     = next+1;
            k        = k+2;         /*  Skip over connected edge */

        else
            k = k+1;                /*  Node has only 1 connection - will be picked up above */

        end
    end
    ncc          = next-1; 
    c1(next:end) = []; 
    c2(next:end) = [];


    /*  Plot the contours
     * If an output is required, extra sorting of the
     * contours is necessary for CLABEL to work.    */
    if nargout > 0 && false

        /*  Form connectivity for the contour, connecting 
         * its edges (rows in cc) with its vertices. */
        ndx = repmat(1,nume,1);
        n2e = 0*penew;
        for k = 1:ncc
            /*  Vertices */
            n1 = c1(k); n2 = c2(k);
            /*  Connectivity */
            n2e(n1,ndx(n1)) = k; ndx(n1) = ndx(n1)+1;
            n2e(n2,ndx(n2)) = k; ndx(n2) = ndx(n2)+1;
        end
        bndn = n2e(:,2)==0;         /*  Boundary nodes */

        bnde = bndn(c1)|bndn(c2);   /*  Boundary edges */


        /*  Alloc some space */
        tmpv = repmat(0,1,ncc);

        /*  Loop through the points at the current contour level (lev(v))
         * Try to assemble the CS data structure introduced in "contours.m"
         * so that clabel will work. Assemble CS by "walking" around each 
         * subcontour segment contiguously. */
        ce    = 1;
        start = ce;
        next  = 2;
        cn    = c2(1);
        flag  = false(ncc,1);        
        x     = tmpv; x(1) = penew(c1(ce),1);
        y     = tmpv; y(1) = penew(c1(ce),2);
        for k = 1:ncc

            /*  Checked this edge */
            flag(ce) = true;

            /*  Add vertices to patch data */
            x(next) = penew(cn,1);
            y(next) = penew(cn,2);
            next    = next+1;

            /*  Find edge (that is not ce) joined to cn */
            if ce==n2e(cn,1)
                ce = n2e(cn,2);
            else
                ce = n2e(cn,1);
            end

            /*  Check the new edge */
            if (ce==0)||(ce==start)||(flag(ce))     

                /*  Plot current subcontour as a patch and save handles */
                x   = x(1:next-1);
                y   = y(1:next-1);
                z   = repmat(lev(v),1,next);
                h   = [h; patch(" Xdata ",[x,NaN]," Ydata ",[y,NaN]," Zdata ",z, ...
                                " Cdata ",z," facecolor "," none "," edgecolor "," flat ")]; hold on      

                /*  Update the CS data structure as per "contours.m"
                 * so that clabel works */
                c = horzcat(c,[lev(v), x; next-1, y]);

                if all(flag)    /*  No more points at lev(v) */

                    break
                else            /*  More points, but need to start a new subcontour */


                    /*  Find the unflagged edges */
                    edges = find(~flag);
                    ce    = edges(1);
                    /*  Try to select a boundary edge so that we are 
                     * not repeatedly running into the boundary */
                    for i = 1:length(edges)
                        if bnde(edges(i))
                            ce = edges(i); break
                        end
                    end
                    /*  Reset counters */
                    start = ce;
                    next  = 2;
                    /*  Get the non bnd node in ce */
                    if bndn(c2(ce))
                        cn = c1(ce);
                        /*  New patch vectors */
                        x = tmpv; x(1) = penew(c2(ce),1);
                        y = tmpv; y(1) = penew(c2(ce),2);
                    else
                        cn = c2(ce);
                        /*  New patch vectors */
                        x = tmpv; x(1) = penew(c1(ce),1);
                        y = tmpv; y(1) = penew(c1(ce),2);
                    end                    

                end

            else                            
                /*  Find node (that is not cn) in ce */
                if cn==c1(ce)
                    cn = c2(ce);
                else
                    cn = c1(ce);
                end
            end

        end

    else        /*  Just plot the contours as is, this is faster... */


        z = repmat(lev(v),2,ncc)+.01*abs(lev(v));

        h = [h; patch(" Xdata ",[penew(c1,1),penew(c2,1)]^t, ...
              " Ydata ",[penew(c1,2),penew(c2,2)]^t, ...
              " Zdata ",z," Cdata ",z," Parent ",ax,varargin[:])]; 

        hold on

    end

end

/*  Assign outputs if needed */
if nargout>0
    cout = c;
    hout = h;
end

return


/* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
}

function [e , eINt , e2t ] = tricontour>mkcon(p,double t) {

numt = size(t,1);
vect = 1:numt;

/*  DETERMINE UNIQUE EDGES IN MESH */

e       = [t(:,[1,2]); t(:,[2,3]); t(:,[3,1])];             /*  Edges - not unique */

vec     = (1:size(e,1))^t;                                   /*  List of edge numbers */

[e,j,j] = unique(sort(e,2)," rows ");                         /*  Unique edges */

vec     = vec(j);                                           /*  Unique edge numbers */

eINt    = [vec(vect), vec(vect+numt), vec(vect+2*numt)];    /*  Unique edges in each triangle */


/*  DETERMINE EDGE TO TRIANGLE CONNECTIVITY */

/*  Each row has two entries corresponding to the triangle numbers
 * associated with each edge. Boundary edges have one entry = 0. */
nume = size(e,1);
e2t  = repmat(0,nume,2);
ndx  = repmat(1,nume,1);
for k = 1:numt
    /*  Edge in kth triangle */
    e1 = eINt(k,1); e2 = eINt(k,2); e3 = eINt(k,3);
    /*  Edge 1 */
    e2t(e1,ndx(e1)) = k; ndx(e1) = ndx(e1)+1;
    /*  Edge 2 */
    e2t(e2,ndx(e2)) = k; ndx(e2) = ndx(e2)+1;
    /*  Edge 3 */
    e2t(e3,ndx(e3)) = k; ndx(e3) = ndx(e3)+1;
end

return

}