function h = pdeplot3D( varargin )
%PDEPLOT3D Plot mesh and solution for 3-D PDE
%   H=PDEPLOT3D(PDEM) plots the mesh for a 3-D PDE defined by the PDEModel
%   object PDEM. The handles to the plotted graphics objects are
%   returned in the optional output argument H.
%
%   H = PDEPLOT(MESH) plots the mesh defined as a FEMesh object. MESH can 
%   be a property of a 3-D model object of type PDEModel or femodel, or a 
%   property of a 3-D results object.
%
%   H=PDEPLOT3D(NODES, ELEMENTS) plots the mesh for a 2-D PDE defined by the 
%   NODES and ELEMENTS matrices. NODES is a 3-by-NumNodes matrix representing 
%   nodal coordinates. ELEMENTS is a 4-by-NumElements or 10-by-NumElements 
%   matrix representing element connectivity in terms of node IDs. 
%
%   H=PDEPLOT3D(...,'Name1',Value1, 'Name2',Value2,...) plots the solution
%   to a 3-D PDE over the specified mesh or elements. 
%
%   The Name-Value pairs are used to specify the solution values and
%   plotting styles are as follows:
%
%   Name        Value/{Default}   Description
%   ----------------------------------------------------------------
%   ColormapData            data          - Colormap plot.
%                                           ColormapData is a vector of nodal 
%                                           solution values.
%   Mesh                    {off} | on    - Display mesh edges. This option
%                                           is used in conjunction with
%                                           ColormapData. It is ignored if
%                                           ColormapData it not specified.
%   NodeLabels              {off} | on    - Display node labels on the 
%                                           mesh boundary. Use in
%                                           conjunction with FaceAlpha.
%   ElementLabels           {off} | on    - Display element labels on the 
%                                           mesh boundary. Use in
%                                           conjunction with FaceAlpha.
%   FaceAlpha               data          - Transparency of faces.
%                                           Scalar in the range [0, 1],
%                                           default = 1
%   Deformation             Displacement  - For a structural analysis model, 
%                                           use this option to plot the
%                                           deformed shape. Value must be
%                                           the displacement structure,
%                                           available in
%                                           StaticStructuralResults object.
%  DeformationScaleFactor  data           - Defines the scaling factor for
%                                           plotting deformed shape. Use
%                                           this option in conjunction with
%                                           Deformation to override the
%                                           default. Default value is
%                                           computed internally based on
%                                           the geometry dimension and
%                                           magnitude of deformation.
% 
%
%   Example 1: Plot a 3-D mesh.
%     pdem = createpde();
%     importGeometry(pdem,'Block.stl');
%     generateMesh(pdem,'Hmax',20);
%     pdeplot3D(pdem,'NodeLabels','on', 'FaceAlpha',0.5);
%
%   Example 2: Solve a cantilever beam problem and plot the deformed shape.
%     pdem = createpde('structural','static-solid');
%     importGeometry(pdem,'SquareBeam.stl');
%     pdegplot(pdem,'FaceLabels','on')
%     generateMesh(pdem);
%     structuralProperties(pdem,'PoissonsRatio',0.3,'YoungsModulus',210E3);
%     structuralBC(pdem,'Face',6,'Constraint','fixed');
%     structuralBoundaryLoad(pdem,'Face',5,'SurfaceTraction',[0,0,-2]);
%     R = solve(pdem);
%     % Plot deformed shape with von Mises stress as contours 
%     figure    
%     pdeplot3D(pdem,'ColorMapData',R.VonMisesStress,'Deformation',R.Displacement)
%     % Plot deformed shape with a scale factor of 100 and von Mises Stress as contours
%     figure
%     pdeplot3D(pdem,'ColorMapData',R.VonMisesStress,'Deformation',...
%                   R.Displacement,'DeformationScaleFactor',100)
%
%   See also: CREATEPDE, pde.PDEMODEL, PDEGPLOT, pde.FEMESH

%   Copyright 2014-2024 The MathWorks, Inc.

% get axis if passed in 
[ax,varargin,nargin] = axescheck(varargin{:});
if ~isempty(ax) && ~isa(ax,'matlab.graphics.axis.Axes')
    error(message('pde:pdegplot:BadAxesType')); 
end
isUserAx = ~isempty(ax);

if nargin > 0
    [varargin{:}] = convertStringsToChars(varargin{:});
end

nargs = nargin;
if nargs < 1
  error(message('pde:pdeplot:nargin'))
end

thepde = [];
nodeAndElementDef = false;
if isa(varargin{1}, 'pde.EquationModel')
  if rem(nargs,2)~=1
        error(message('pde:pdeplot:NoParamPairs'))
  end
  thepde = varargin{1};
  if ~isa(thepde.Mesh, 'pde.FEMesh')
      error(message('pde:pdeplot:UnmeshedPdeModel'))
  elseif thepde.IsTwoD
      error(message('pde:pdeplot:NotThreeD'))
  else
      themsh = thepde.Mesh;
      [p,~,t] = themsh.meshToPet();    
      varargin = {varargin{2:end}};      
  end
elseif isa(varargin{1}, 'pde.FEMesh') 
  themsh = varargin{1};
  if rem(nargs,2)~=1
      error(message('pde:pdeplot:NoParamPairs'))
  end
  if size(themsh.Nodes,1) ~= 3
      error(message('pde:pdeplot:NotThreeD'))
  else      
      [p,~,t] = themsh.meshToPet();    
      varargin = {varargin{2:end}};      
  end  
elseif ~isnumeric(varargin{1}) % something other than p,e,t
    error(message('pde:pdeplot:InvalidArgs'))
else
    if nargs < 2
        error(message('pde:pdeplot:InvalidArgs'))
    elseif rem(nargs,2)~=0
        error(message('pde:pdeplot:NoParamPairs'))
    end
    p = varargin{1};
    t = varargin{2};
    nodeAndElementDef = true;
    validateattributes(varargin{1},{'numeric'},{'real', 'finite', 'nonsparse', 'nonnan'},'pdeplot3D',inputname(1));
    validateattributes(varargin{2},{'numeric'},{'real', 'integer'},'pdeplot3D',inputname(2));
    varargin = {varargin{3:end}};    
end

if size(p,1) ~= 3
   error(message('pde:pdeplot:NotThreeD'))
end   

parser = inputParser;
addParameter(parser,'colormapdata', [], @isValidU);
addParameter(parser,'flowdata', [], @isnumeric);
addParameter(parser,'Deformation', []);
addParameter(parser,'DeformationScaleFactor', [], @isnumeric);
addParameter(parser,'NodeLabels', 'off', @isValidNdLabelOption);
addParameter(parser,'ElementLabels', 'off', @isValidElLabelOption);
addParameter(parser,'Mesh', 'off', @isValidMeshDispOption);
addParameter(parser,'FaceAlpha', 1, @isnumeric);
addParameter(parser,'EdgeColor', 'k');
addParameter(parser,'FaceColor', 'cyan');
parse(parser,varargin{:});
colormapdata = parser.Results.colormapdata;
flowdata = parser.Results.flowdata;
faceAlpha = parser.Results.FaceAlpha;
showMesh = parser.Results.Mesh;
showNodeLabels = parser.Results.NodeLabels;
showElemLabels = parser.Results.ElementLabels;
deformation = parser.Results.Deformation;
scaleFactor = parser.Results.DeformationScaleFactor;
faceColor = parser.Results.FaceColor;
edgeColor = parser.Results.EdgeColor;

if ~isnumeric(faceAlpha) || ~isscalar(faceAlpha) || ~isfinite(faceAlpha) || faceAlpha < 0 || faceAlpha > 1
    error(message('pde:pdeplot:InvalidAlphaOption'))
end

if(strcmpi(showElemLabels, 'on') && nodeAndElementDef && size(p,1) == 3) 
    error(message('pde:pdeplot:ElemLabelsNeedsCompleteMesh'));
end

if(size(t,1) == 5 || size(t,1) == 11)
    t(end,:)=[];
end   

numElemNodes = size(t,1);
if(numElemNodes ~= 4 && numElemNodes ~= 10)
  error(message('pde:pdeModel:invalidT'));
end
numNodes = size(p,2);
if ~isempty(colormapdata)
  if(length(colormapdata) ~= numNodes)
    error(message('pde:pdeplot:colormapdataLength'));
  end
end

if isempty(ax)
  ha = newplot;
else
  ha = ax;
  if ~ishold(ha)
    cla(ha,'reset');
  end
end
hf = ancestor(ha,'figure');

theme = updateTheme(hf);
set(ha,'ClippingStyle','rectangle');
hold(ha,'on');

bbox = [min(p(1,:)) max(p(1,:));
        min(p(2,:)) max(p(2,:));
        min(p(3,:)) max(p(3,:))];
    
if ~isempty(deformation) && isValidDeformation(deformation)
    if isempty(scaleFactor)
        MaxDeformationMag = max(sqrt(deformation.ux.^2+deformation.uy.^2+deformation.uz.^2));
        if(MaxDeformationMag > eps())
            scaleFactor = min(bbox(:,2) - bbox(:,1)) /MaxDeformationMag; % Based on lowest bounding box dimension
        else
            scaleFactor = 0;
        end
    end
    p(1,:) = p(1,:) +  scaleFactor*deformation.ux';
    p(2,:) = p(2,:) +  scaleFactor*deformation.uy';
    p(3,:) = p(3,:) +  scaleFactor*deformation.uz';
end

bbox = [min(p(1,:)) max(p(1,:));
        min(p(2,:)) max(p(2,:));
        min(p(3,:)) max(p(3,:))];
 
plotAxisTriad(bbox,ha);


if ~isempty(colormapdata) 
  [ltri] = tetBoundaryFacets(p,t);  
  if(numElemNodes == 10)
      ltri = splitQuadraticTri(ltri);    
  end
  h1=colorbar(ha); ht = hgtransform(ha); colormap(ha,'jet');
  if min(colormapdata) ~= max(colormapdata)
      caxis(ha, [min(colormapdata) max(colormapdata)]);
  end
  
  h2=patch('Faces',ltri, 'Vertices', p', 'FaceVertexCData', colormapdata(:), ...
    'AmbientStrength', .75,  ...
    'EdgeColor', 'none', 'FaceColor', 'interp', 'parent',ht,'FaceAlpha',faceAlpha, 'Clipping','off');
  if strcmpi(showMesh, 'on')   
    [ltri, lp] = tetBoundaryFacets(p,t(1:4,:));  
    tr = triangulation(ltri, lp);
    tre = (tr.edges())';
    x = lp(:,1);
    y = lp(:,2);
    z = lp(:,3);
    xt=x(tre); xt(3,:) = NaN; 
    yt=y(tre); yt(3,:) = NaN;
    zt=z(tre); zt(3,:) = NaN;
    plot3(ha,xt(:),yt(:),zt(:),'-k');    
  end
  if(strcmpi(showNodeLabels, 'on'))
    plotNodeLabels(p, t, ha);
  end
  if(strcmpi(showElemLabels, 'on') && ~isempty(thepde))      
      plotElementLabels(thepde.Mesh, ha);
  end     

  if nargout==1
    h = [h1 h2];
  end
elseif ~isempty(flowdata)
    if(length(flowdata(:)) ~= 3*size(p,2))
       error('Length of flowdata must be 3*number of points');
    end
    if isvector(flowdata)
      flowdata = reshape(flowdata, numNodes, 3);
    end
    quiver3(ha,p(1,:)', p(2,:)', p(3,:)', flowdata(:,1), flowdata(:,2), flowdata(:,3));
else
  %Plot the boundary triangulation:
  [ltri, lp] = tetBoundaryFacets(p,t(1:4,:));  
  hf = ancestor(ha,'figure'); % to set renderer, must get figure, but that may not be direct parent
  hh=trisurf(ltri, lp(:,1),lp(:,2),lp(:,3), ...
    'FaceColor',faceColor, 'EdgeColor',edgeColor,'FaceAlpha', faceAlpha, 'Clipping','off','Parent',ha);
  if nargout==1
    h = hh;
  end    
  if(strcmpi(showNodeLabels, 'on'))
    plotNodeLabels(p,t,ha);
  end
  if(strcmpi(showElemLabels, 'on'))
     plotElementLabels(themsh,ha);
  end     
end
axis(ha,'equal','tight','off');
view(ha,30,30);
hold(ha,'off');
end


function plotNodeLabels(p,t,ha)
    tri = tetBoundaryFacets(p,t);      
    vispts = unique(tri(:));   
    warnState = warning('off', 'MATLAB:triangulation:PtsNotInTriWarnId');    
    tr = triangulation(tri(:,1:3),p');
    if size(tri,2) >= 6
        tr2 = triangulation(tri(:,4:6),p');
    end
    warning(warnState);
    [~, ric] = tr.incenter();  
    ric = mean(ric);
    vn = tr.vertexNormal(vispts);
    if size(tri,2) >= 6
        vn = vn + tr2.vertexNormal(vispts);
    end
    vn = (vn.*(0.15*ric))';
    px = p(1,vispts)+ vn(1,:);
    py = p(2,vispts)+ vn(2,:);
    pz = p(3,vispts)+ vn(3,:);
    labels = 'n' + string(1:max(vispts)); labels = cellstr(labels');
    labels = labels(vispts);
    pos = single([px ; py; pz]);
    % set label color
    theme = updateTheme(ha.Parent);
    if strcmpi(theme.BaseColorStyle,'dark')
        labColor = uint8([255 255 255 255]');
    else
        labColor = uint8([0 0 0 255]');
    end
    nt = matlab.graphics.primitive.world.Text('VertexData', pos, 'String', labels,...
        'HorizontalAlignment','center', 'Clipping','on','Interpreter','none', ...
        'Parent', ha, 'ColorData', labColor);
    nt.Font.Name = 'Helvetica';
    nt.Font.Size=8;
end

function plotElementLabels(msh,ha)
    ma = msh.MeshAssociation;
    fa = ma.FaceAssociativity;
    facets = cell2mat(fa');
    n = msh.Nodes';
    e = msh.Elements(1:4,:)';
    lfc = ma.FaceCodes(:,1:3); % Linear Face codes
    facets = facets';
    i = sub2ind(size(e), facets(:,1), lfc(facets(:,2),1));
    j = sub2ind(size(e), facets(:,1), lfc(facets(:,2),2));
    k = sub2ind(size(e), facets(:,1), lfc(facets(:,2),3));
    ei = e(i);
    ej = e(j);
    ek = e(k);
    nf = [ei(:), ej(:), ek(:)];
    warnState = warning('off', 'MATLAB:triangulation:PtsNotInTriWarnId');    
    tr = triangulation(nf,n);
    warning(warnState);
    [ic, ric] = tr.incenter();
    fn = faceNormal(tr);
    ic = ic + (fn.*ric)*0.15;
    labels = 'e' + string(1:max(facets(:,1))); labels = cellstr(labels');
    labels = labels(facets(:,1));
    pos = single(ic');
    % set label color
    theme = updateTheme(ha.Parent);
    if strcmpi(theme.BaseColorStyle,'dark')
        labColor = uint8([255 255 255 255]');
    else
        labColor = uint8([0 0 0 255]');
    end
    et = matlab.graphics.primitive.world.Text('VertexData', pos, 'String', labels,...
        'HorizontalAlignment','center', 'Clipping','on','Interpreter','none', ...
        'Parent', ha, 'ColorData', labColor);
    et.Font.Name = 'Helvetica';
    et.Font.Size=8;
end


function ok = isValidU(udata)
    if ~isreal(udata) || issparse(udata) || any(~isfinite(udata))
        error(message('pde:pdeplot:InvalidColormapdata'));   
    end     
    ok = true;
end

function ok = isValidMeshDispOption(dispopt)
    if ~ischar(dispopt)
       error(message('pde:pdeplot:GenericMustBeOffOn','mesh'));   
    end
    validatestring(dispopt,{'on', 'off'}) ;
    ok = true;
end

function ok = isValidNdLabelOption(labelopt)
    if ~ischar(labelopt)
       error(message('pde:pdeplot:GenericMustBeOffOn', 'NodeLabels'));   
    end
    validatestring(labelopt,{'on', 'off'}) ;
    ok = true;
end

function ok = isValidElLabelOption(labelopt)
    if ~ischar(labelopt)
       error(message('pde:pdeplot:GenericMustBeOffOn', 'ElementLabels'));   
    end
    validatestring(labelopt,{'on', 'off'}) ;
    ok = true;
end

function ok = isValidDeformation(deformation)
    if (isstruct(deformation)||isa(deformation,"FEStruct")) ...
        && isscalar(deformation) ...
        && isFieldProp(deformation,"ux") && isFieldProp(deformation,"uy") ...
        && isFieldProp(deformation,"uz")
        ok = true;
        return
    end
    ok = false;
end

function ok = isFieldProp(festruct,field)
    ok = isfield(festruct,field) || isprop(festruct,field);
end

function t4=splitQuadraticTri(t)
t4Nodes = [1 4 6; 4 5 6; 4 2 5; 6 5 3];
t4 = [t(:,t4Nodes(1,:)); t(:,t4Nodes(2,:)); t(:,t4Nodes(3,:)); t(:,t4Nodes(4,:))];
end

function theme = updateTheme(fg)
% Respond to theme changes
if isempty(fg.Theme)
    % only run if figure theme is empty
    matlab.graphics.internal.themes.figureUseDesktopTheme(fg);
end
% if figure has a theme, return it, else return light
theme = matlab.graphics.internal.themes.lightTheme; 
if ~isempty(fg.Theme)
    theme = fg.Theme;
end
end