function mesh = findfixednodes( mesh )
%mesh = findfixednodes( mesh )    Find three nodes whose degrees of freedom
%   can be contrained without distorting or rotating the mesh too much.

%First find two nodes where the direction from one to the other is nearly
%parallel to the x axis.

    method = 0;
    
    if method==0
        mesh.globalProps.fixedDFs = [];
        return;
    end
    
    numnodes = size( mesh.nodes, 1 );
    
    switch method
        case 1
            [minx,ni] = min( mesh.nodes(:,1) );
            [maxx,nj] = max( mesh.nodes(:,1) );
            v = mesh.nodes(nj,:) - mesh.nodes(ni,:);
            minangle = atan2( v(2), v(1) );
            minabsangle = abs(minangle);
        case 2
            minabsangle = pi*3;
            minangle = minabsangle;
            for i=1:numnodes-1
                for j=i+1:numnodes
                    v = mesh.nodes(i,:) - mesh.nodes(j,:);
                    angle = atan2( v(2), v(1) );
                    absangle = abs(angle);
                    if absangle < minabsangle
                        minabsangle = absangle;
                        minangle = angle;
                        ni = i;
                        nj = j;
                    end
                end
            end
        case 3
            numedges = size( mesh.edgeends, 1 );
            edgedirs = zeros( 1, numedges );
            for i=1:numedges
                v = mesh.nodes(mesh.edgeends(i,2),:) - mesh.nodes(mesh.edgeends(i,1),:);
                edgedirs(i) = atan2( v(2), v(1) );
            end
          % [edgesort,i] = sort(edgedirs,2);
            [minabsangle,ei] = min( abs( edgedirs ) );
            minangle = edgedirs(ei);
            ni = mesh.edgeends(ei,1);
            nj = mesh.edgeends(ei,2);
        otherwise
    end
    if minabsangle > 0.000001
        fprintf( 1, 'Method %d, min angle %f for nodes %d and %d.\n', ...
            method, minangle, ni, nj );
        mesh = rotatemesh( mesh, -minangle, 'Z' );
        mesh.globalProps.prismnodesvalid = false;
        mesh = makeprismsvalid( mesh );
        v = mesh.nodes(nj,:) - mesh.nodes(ni,:);
        newminangle = atan2( v(2), v(1) );
        fprintf( 1, 'New min angle %f [%f,%f,%f] [%f,%f,%f].\n', ...
            newminangle, mesh.nodes(ni,:), mesh.nodes(nj,:) );
    end

    fni2 = ni*2;
    fni1 = fni2-1;
    fnj = nj*2;
    dfi1base = fni1*3;
    dfi2base = fni2*3;
    dfjbase = fnj*3;
    mesh.globalProps.fixedDFs = ...
        [ dfi1base-2, dfi1base-1, dfi1base, dfi2base-2, dfi2base-1, dfjbase-1 ];
end
