function out = pdesmech(p,t,c,u,NVArgs)
%PDESMECH  Calculates structural mechanics tensor functions.
%     PDESMECH is not recommended. Use the PDE app instead. 
%
%     OUT=PDESMECH(P,T,C,U,P1,V1,...) computes tensor expressions such
%     as stresses and strains for plane stress and plane strain structural
%     mechanics applications.
%     P and T are the mesh point and triangle matrices, respectively,
%     C is the PDE equation's C coefficient, U = PDE solution [u; v] of
%     length 2*number of nodes.
% 
%     Valid Property/value pairs include
% 
%     tensor          'ux' | 'uy' | 'vx' | 'vy' | 'exx' | 'eyy' | 'exy' |
%                     'sxx' | 'syy' | 'sxy' | 'e1' | 'e2' | 's1' | 's2' |
%                     {'vonmises'}
%     application     {'ps'} | 'pn'
%     nu              {0.3} numerical value or string expression

%     Copyright 1994-2023 The MathWorks, Inc.

arguments
    p
    t
    c
    u
    NVArgs.tensor {mustBeNonempty, mustBeMember(NVArgs.tensor, ...
        ["ux", "uy", "vx", "vy", "exx", "eyy", "exy", "sxx", "syy", "sxy", ...
        "e1", "e2", "s1", "s2", "vonmises", "von Mises"])} = "vonmises"
    NVArgs.application {mustBeNonempty, mustBeMember(NVArgs.application, ...
        ["ps", "pn"])} = "ps"
    NVArgs.nu (1,:) {mustBeA(NVArgs.nu, ["double", "single", "char", "string"])} = 0.3
end

tensor      = NVArgs.tensor;
application = NVArgs.application;
nu          = NVArgs.nu;
if ischar(nu) || isstring(nu)
    try
        nu = pdetxpd(p,t,u,nu);
    catch
        nu = NaN;
    end
end
if nu<0 || nu>0.5 || ~isreal(nu) || isnan(nu)
    error(message('pde:pdesmech:InvalidNu'))
end

if length(u)~=2*size(p,2)
    error(message('pde:pdesmech:LengthU'))
end

[ux,uy] = pdegrad(p,t,u);
[u1,u2] = pdeeigx(ux(1,:),(ux(2,:) + uy(1,:))/2, ...
    (ux(2,:) + uy(2,:))/2,uy(2,:));
[sx,sy] = pdecgrad(p,t,c,u);
[s1,s2] = pdeeigx(sx(1,:),(sx(2,:) + sy(1,:))/2, ...
    (sx(2,:) + sy(1,:))/2,sy(2,:));

switch tensor
    case "ux"
        out = ux(1,:);
    case "uy"
        out = uy(1,:);
    case "vx"
        out = ux(2,:);
    case "vy"
        out = uy(2,:);
    case "exx"
        out = ux(1,:);
    case "eyy"
        out = uy(2,:);
    case "exy"
        out = ux(2,:) + uy(1,:);
    case "sxx"
        out = sx(1,:);
    case "syy"
        out = sy(2,:);
    case "sxy"
        out = (sx(2,:) + sy(1,:))/2;
    case "e1"
        out = u1';
    case "e2"
        out = u2';
    case "s1"
        out = s1';
    case "s2"
        out = s2';
    case {"vonmises", "von Mises"}
        if strcmp(application, "ps")
            out = sqrt(s1.*s1 + s2.*s2 - s1.*s2)';
        elseif strcmp(application, "pn")
            out = sqrt((s1.*s1 + s2.*s2).*(nu.*nu - nu + 1)' + ...
                s1.*s2.*(2*nu.*nu - 2*nu - 1)')';
        end
end
end