function[L,Astr,Estr,bodystr,Fbar,DBCNodeCode,DBCNodeValue,...
         nel,nen,nnp,eltype,nip,xhat,nelpp,...
         xexact,uexact] = me180_lab06_input_deck()

    % Length of the domain.
    L=10;

    % The following are strings that can be used by the 'eval' command to
    % calculate the area, stiffness, and body forces.
    % I will always evaluate them at a point 'xe'.
    Astr = 'pi*(0.1^2)';
    Estr = '7e10*(1+xe)';
    bodystr = '500 + 10*xe';

    % Neumann BCs.
    Fbar = [0; 500];

    % For Dirichelt BCs, we can't just have a list of values; we need to know which
    % nodes (or positions) have Dirichlet BCs, and then for each node (or position),
    % we need to know the value of the Dirichlet BC.
    %
    % Since we don't know which node corresponds to what physical location,
    % the positions in the DBC structures correspond to the left and right of the
    % domain.
    DBCNodeCode = [1; 0];
    DBCNodeValue = [0; 0];

    % number of elements
    nel=2;

    % element type: 1=linear
    %               2=quadratic
    eltype=2;

    % Calculations of some other quantities based on the element type;
    % these are not to be modified by the user.
    % nen = number of element nodes (i.e., the number of nodes
    %                                      per element)
    % nnp = number of nodal points (in the whole domain)
    if (eltype == 1)
        nen = 2;
        nnp = nel+1;
    elseif (eltype == 2)
        nen = 3;
        nnp = (nel*2)+1;
    end

    % number of integration points
    nip = 3;

    % mesh type: 0=typical
    %            1=random
    meshtype = 1;

    % Calculate the vector of nodal positions.
    xhat = Buildxhat(meshtype,L,nnp);

    % Number of element plot points: this is used to ensure that the element
    % curvature is plotted appropriately.
    nelpp = 10;

    % Create the exact solution
    A = pi*(0.1^2);
    Ebar = 7e10;
    b0 = 500;
    b1 = 10;
    x = linspace(0,L);
    xexact = x;
    term1 = -log(1+x) + (1+x);
    term2 = (x.^2).*log(1+x) ...
           -((1+x).^2).*log(1+x) ...
           +((1+x).^2)/2 ...
           +2*(1+x).*log(1+x) ...
           -2*(1+x);
    term3 = log(1+x);
    c1 = 500 + b0*L + b1/2*L^2;
    c2 = 1/(A*Ebar)*(b0 - 3/4*b1);
    uexact = 1/(A*Ebar)*(-b0*term1 - b1/2*term2 + c1*term3) + c2;


end

function[xhat] = Buildxhat(type,L,nnp)

    if (type == 0)
        xhat = linspace(0,L,nnp);
    else
        used = zeros(nnp,1);
        nums = round(nnp.*rand(nnp*100,1));
        i = 1;
        j = 0;
        h = L/(nnp-1);
        while (1)
            nodenum = nums(i);
            if (nodenum == 0)
                nodenum = 1;
            end
            i = i+1;
            if (used(nodenum)== 1)
                continue;
            end
            used(nodenum) = 1;
            xhat(nodenum) = j*h;
            j = j+1;
            if (j == nnp)
                break;
            end
        end
    end
end


