classdef Shape2D < shape.Shape
    % SHAPE2D is the  superclass for 2d shapes
    %
    % See also <a href="matlab:help shape.Shape2D.Rectangle">rectangle</a>, <a href="matlab:help shape.Shape2D.Circle">circle</a>

    % Copyright 2022-2023 The MathWorks, Inc.
    properties (Hidden,SetAccess=protected)
        FillPolygons
        HolePolygons
        FeaturePolygons
    end

    properties(Hidden,Dependent)
        ShapeVertices
    end

    properties(Hidden)
        InternalPolyShape
        SliverTol = 1e-12;
    end

    methods
        function obj = Shape2D(vertices)
            validateattributes(vertices,{'numeric'},{'nonempty','2d',...
                'real'});
            if size(vertices,2)>3 ||isscalar(vertices)
                error(message('antenna:shapeerrors:InvalidNumCols'));
            end
            if any(any(isinf(vertices)))
                error(message('antenna:shapeerrors:InvalidInpInf'));
            end
            obj.InternalPolyShape   =polyshape(vertices(:,1:2),'Simplify',false);
            obj.FillPolygons    = {};
            obj.HolePolygons    = {};
            obj.FeaturePolygons = [];
            obj.ShapeDimension = '2D';
        end
    end

    methods

        function propVal = get.ShapeVertices(obj)
            % Check if Shape is still on XY
            if isShapeOnXY(obj)
                propVal     = obj.InternalPolyShape.Vertices;
                if size(propVal,2)==2
                    propVal(:,3) = 0;
                end
            else
                G = getGeometry(obj);
                propVal = G.BorderVertices;
            end
        end

        function propVal = get.FillPolygons(obj)
            propVal = obj.FillPolygons;
        end

        function propVal = get.HolePolygons(obj)
            propVal = obj.HolePolygons;
        end
        

        % Visualization Methods

        function varargout = plot(obj,varargin)
            %PLOT   Plot the boundary of the shape
            %
            % P = plot(obj,varargin) plots the boundaries of the shape and
            % also returns the line handle.
            %
            % Example
            %
            % % Example1:Create a rectangular shape and plot it.
            %
            % r  = shape.Rectangle;
            % p  = plot(r);
            %
            % See also <a href="matlab:help shape.Shape2D.show">show</a>, <a href="matlab:help shape.Shape2D.mesh">mesh</a>


            if isempty(obj)
                error(message("antenna:shapeerrors:EmptyShape"));
            end

            if getHasStructureChanged(obj)
                createGeometry(obj);
            end
            [x,y] = boundary(obj.InternalPolyShape);
            % Plot final structure
            if nargout==0
                plot(x,y,varargin{:});
            else
                varargout{1} = plot(x,y,varargin{:});
            end

            axis equal
        end


        function varargout = extrudeLinear(shapein,height,varargin)
            % Extrudes the given 2d shape linearly along the given direction for a 
            % given height. Default extrusion is along z-axis.
            %
            % N-V Pairs:
            % Twist       - Twists the given shape by the specified angle for each segment over the height of
            %               extrusion
            % Scale       - For a scalar input,shape is scaled by given value over
            %               height of extrusion. For a vector input, shape is scaled
            %               in x-direction by first element of vector and shape is
            %               scaled in y-direction by second element of vector over the
            %               height of extrusion
            % NumSegments - Number of divisions over the height of extrusion
            % Caps        - Logical value determines start and end points of extrusion to be
            %               closed or opened. For logical true,caps are
            %               present and for logical false caps are removed.
            % Direction   - Direction in which linear extrusion is done.
            %
            % % Examples:
            %
            % % Example 1:
            % % Extrude a rectangle linearly for a height of 1 m.
            %
            % rect = shape.Rectangle;
            % extrudeLinear(rect,1);
            %
            % % Example 2:
            % % Extrude a circle by applying scaling along height of
            % % extrusion
            %
            % circ = shape.Circle;
            % extrudeLinear(circ,1,Scale = 2);
            %
            % % Example 3:
            % % Extrude a rectangle by increasing the segments along height
            % % and by applying a twisting angle of 1 degree.
            %
            % rect = shape.Rectangle;
            % extrudeLinear(rect,1,Twist = 1,NumSegments = 30);

  
            % Get input parameters
            parserObj=inputParser;
            addParameter(parserObj,'Twist',0);
            addParameter(parserObj,'Direction',[0 0 1]);
            addParameter(parserObj,'Scale',1);
            addParameter(parserObj,'NumSegments',1);
            addParameter(parserObj,'Caps',false);
            addParameter(parserObj,'Parent',[]);
            parse(parserObj,varargin{:});
            Twist   = parserObj.Results.Twist;
            Direction   = parserObj.Results.Direction;
            Scale   = parserObj.Results.Scale;
            NumSegments   = parserObj.Results.NumSegments;
            Caps   = parserObj.Results.Caps;
            ParentObj   = parserObj.Results.Parent;


            % Validate inputs
            nameOfFunction = 'extrudeLinear';
            validateattributes(shapein,{'shape.Shape2D'},{'nonempty','scalar'},      ...
                nameOfFunction,'Shape',1);

            validateattributes(height,{'numeric'},{'nonempty','scalar',      ...
                'nonnan','finite','real','positive'},   ...
                nameOfFunction,'Height',2);

            validateattributes(Twist,{'numeric'},{'nonempty','scalar',      ...
                'nonnan','finite','real','nonnegative','>=',0,'<=',360},   ...
                nameOfFunction,'Twist');

            validateattributes(Direction,{'numeric'},{'nonempty','vector',      ...
                'nonnan','finite','real','numel',3},   ...
                nameOfFunction,'Direction');

            if isscalar(Scale)
                validateattributes(Scale,{'numeric'},{'nonempty','scalar',      ...
                    'nonnan','finite','real','positive'},   ...
                    nameOfFunction,'Scale');
            else
                validateattributes(Scale,{'numeric'},{'nonempty','vector',      ...
                    'nonnan','finite','real','positive','numel',2},   ...
                    nameOfFunction,'Scale');
            end

            validateattributes(NumSegments,{'numeric'},{'nonempty','scalar',      ...
                'nonnan','finite','real','integer','positive'},   ...
                nameOfFunction,'NumSegments');
            validateattributes(Caps,{'logical'},{'scalar','binary'}, ...
                nameOfFunction,'Caps');

            if ~isempty(ParentObj)
                validateattributes(ParentObj,["matlab.ui.Figure","matlab.graphics.axis.Axes","matlab.ui.control.UIAxes"],{'scalar'},...
                    nameOfFunction,'Parent');
            end
            checkParentInput(shapein,ParentObj);
            
            % Get vertices of 2d shape
            vert = getShapeVertices(shapein);

            % Error out if 2d shape is not on x-y plane
            if any(vert(:,3)>0)
                     error(message('antenna:shapeerrors:Unsupported', ...
                    'ExtrudeLinear','shapes not on XY-plane '));  
            end

            % Error out if direction takes unacceptable value
            if all(Direction == 0)
                error(message('antenna:shapeerrors:InvalidDirection'));
            end

            % get multiplier
            mult1=em.internal.authoring.getMultiplier(shapein,vert);
            if isa(shapein,'shape.Rectangle')% Required due to we r calculating exponent rather than number of digits after decimals
                mult2=em.internal.authoring.getMultiplier(shapein,[height,shapein.Center,shapein.Length,shapein.Width],1);
            elseif isa(shapein,'shape.Circle')
                mult2=em.internal.authoring.getMultiplier(shapein,[height,shapein.Center,shapein.Radius],1);
            end
            if ~isa(shapein,'polyshape')&&~isa(shapein,'shape.Polygon')
                if isempty(shapein.Multiplier)
                    mult = max(mult1,mult2);
                else
                    mult=shapein.Multiplier;
                end
            else
                mult=mult1;
            end

            % New Changes
            poly_shape = polyshape(vert(:,1:2));
            % matid = unique(shapein.MultiShapeID);
            % if isempty(matid)
                matid = 1;
            % end
            meshobj = Geometry.linearExtrude(poly_shape,'Height',height,...
               'Steps',NumSegments,'EndRotation',Twist*NumSegments,...
               'EndScale',Scale,'Direction',Direction,'Material',string(matid),'Caps',Caps);
            meshobj = Geometry.meshClean(meshobj);
            [pp,tt] = Geometry.extractMesh(meshobj);
            obj = shape.Custom3D;obj.Points = pp;obj.Triangles=tt(:,1:3)+1;
            % obj.Points = obj.Points./mult;
            obj.Points = round(obj.Points,log10(mult*1e4));
            obj.MultiShapeID = tt(:,4);
            %

            % Show final structure
            if nargout==0
                if isempty(ParentObj)
                    show(obj);
                else
                    show(obj,"Parent",ParentObj);
                end
            else
                varargout{1} = obj;
            end
        end
        
        function varargout = extrudeRotate(shapein,angle,options)
            % extrudeRotate generates extrusion of the given 2d shape by rotating  upto
            % the given angle about Z-Axis
            %
            % N-V Pairs:
            % NumSegments - Number of divisions along the extrusion
            % Caps        - Logical value determines start and end 
            %               extrusion to be closed or opened.
            % Pitch       - Distance along z-axis for one complete rotation
            %
            % % Examples:
            %
            % % Example 1:
            % % Extrude a circle by rotating it by an angle of 90 degree and
            % % by using 30 segments
            %
            % circ = shape.Circle(Center = [2 0]);
            % extrudeRotate(circ,90,NumSegments = 30);
            
            arguments
                shapein (1,1) shape.Shape2D
                angle  {mustBeNumeric,mustBeNonempty,...
                    mustBeScalarOrEmpty,mustBeReal,mustBeNonNan,mustBeFinite...
                    };
                options.NumSegments  {mustBeNumeric,mustBeNonempty,...
                    mustBeScalarOrEmpty,mustBeReal,mustBeNonNan,mustBeFinite,...
                    mustBePositive,mustBeInteger...
                    } =3;
                options.Parent  {mustBeA(options.Parent,...
                    ["matlab.ui.Figure","matlab.graphics.axis.Axes","matlab.ui.control.UIAxes"])}
                options.Caps (1,1) logical = false;
                options.Pitch {mustBeNumeric,mustBeNonempty,...
                    mustBeScalarOrEmpty,mustBeReal,mustBeNonNan,mustBeFinite...
                    }=0;
             end
             
             NumSegments = options.NumSegments;
             Caps = false;
             Pitch = options.Pitch;
             if isfield(options,'Parent')
                 ParentObj = options.Parent;
             else
                 ParentObj =[];
             end
             checkParentInput(shapein,ParentObj);

            
%              
%             NumSegments   = parserObj.Results.NumSegments;
%             ParentObj   = parserObj.Results.Parent;

            if getHasStructureChanged(shapein)
                createGeometry(shapein);
            end
            % get the vertices and generate triangulation
            vert = shapein.ShapeVertices;
            if any(vert(:,3)>0)
                     error(message('antenna:shapeerrors:Unsupported', ...
                    'Extrude','shapes not on XY-plane '));  
            end

            mult1=em.internal.authoring.getMultiplier(shapein,vert);
            if isa(shapein,'shape.Rectangle')% Required due to we r calculating exponent rather than number of digits after decimals
                mult2=em.internal.authoring.getMultiplier(shapein,[angle,shapein.Center,shapein.Length,shapein.Width],1);
            elseif isa(shapein,'shape.Circle')
                mult2=em.internal.authoring.getMultiplier(shapein,[angle,shapein.Center,shapein.Radius],1);
            end
            if ~isa(shapein,'polyshape')&&~isa(shapein,'shape.Polygon')
                mult = max(mult1,mult2);
            else
                mult=mult1;
            end

            vert = vert*mult;
            poly_shape = polyshape(vert(:,1:2));
            matid = unique(shapein.MultiShapeID);
            if isempty(matid)
                matid = 1;
            end
            meshobj = Geometry.circularExtrude(poly_shape,'Angle',angle,...
               'Steps',NumSegments,...
               'Material',string(matid),'Caps',Caps,'Pitch',Pitch);
            meshobj = Geometry.meshClean(meshobj);
            [pp,tt] = Geometry.extractMesh(meshobj);
            obj = shape.Custom3D;obj.Points = pp;obj.Triangles=tt(:,1:3)+1;
            obj.Points = obj.Points./mult;
            obj.Points = round(obj.Points,log10(mult*1e4));
            obj.MultiShapeID = tt(:,4);

            % Show final structure
            if nargout==0
                if isempty(ParentObj)
                    show(obj);
                else
                    show(obj,"Parent",ParentObj);
                end
            else
                varargout{1} = obj;
            end
        end



        function show(obj,varargin)
            %SHOW   Visualize the shape as a filled patch
            %
            % show(obj) plots the shape as a filled region using patches.
            %
            % Example
            %
            % % Example1: Create a circular shape and visualize the filled
            % % regions
            %
            % c  = shape.Circle;
            % show(c);
            %
            % See also <a href="matlab:help shape.Shape2D.plot">plot</a>, <a href="matlab:help shape.Shape2D.mesh">mesh</a>

           [ax,~,colflag]=parseShowInputs(obj,varargin{:});
           checkParentInput(obj,ax);

            if isempty(obj)
                error(message("antenna:shapeerrors:EmptyShape"));
            end

            if getHasStructureChanged(obj)
                createGeometry(obj);
            end
            updateGeometryCache(obj);
            geometry        = getGeometry(obj);
%             antennaColor    = [223, 185, 58]/255;
            if iscell(geometry)
                BorderVertices = cell2mat(cellfun(@(x) x.BorderVertices,...
                    geometry', 'UniformOutput',false));
                [~, mult, unit] = engunits(BorderVertices);
                for m = 1:numel(obj.Geometry)
                    geometry{m}.multiplier = mult;
                    geometry{m}.unit       = unit;
                end
            else
                BorderVertices = geometry.BorderVertices;
                [~, mult, unit] = engunits(BorderVertices);
                geometry.multiplier = mult;
                geometry.unit       = unit;
            end
            ScaledPoints = obj.Points*mult;

            % Unique shape id to know how many shapes are present
            usid=unique(obj.MultiShapeID);
            if ~isempty(usid)
                P(:,1:usid(1))=0;
            end
            trr=cell2mat(geometry.polygons);
            
            [fColor,colid] = getShowParam(obj,colflag);
            if ~isempty(obj.MultiShapeID)
                % This condition is when  boolean operation is performed
                for i=1:size(usid,1)
                    if ~isempty(obj.ColorID)
                        P1=patch(ax,'Vertices',ScaledPoints,'Faces',...
                            trr((obj.MultiShapeID(:,1))==i,1:3),...
                            'FaceColor',colid{i,1},'EdgeColor',...
                            'None','faceAlpha',obj.ColorID{i,2},'EdgeAlpha',obj.ColorID{i,2});
                    else
                         P1=patch(ax,'Vertices',ScaledPoints,'Faces',...
                            cell2mat(geometry.polygons),...
                            'FaceColor','y','EdgeColor',...
                            'None','faceAlpha',1,'EdgeAlpha',1);
                    end
                    P(:,i)=P1;
                    hold on;
                end
                if ~isempty(obj.ColorID)
                    [~]=patch(ax,'Vertices',ScaledPoints,'Faces',cell2mat(geometry.BoundaryEdges),'EdgeColor',obj.ColorID{i,3},'FaceColor','None');
                else
                    [~]=patch(ax,'Vertices',ScaledPoints,'Faces',cell2mat(geometry.BoundaryEdges),'EdgeColor','k','FaceColor','None');
                end
            else
                % When boolean operation is not performed on the shape
                [~] = patch(ax,'Vertices',ScaledPoints,'Faces',...
                    cell2mat(geometry.polygons),'FaceColor',fColor,'EdgeColor',...
                    'None','faceAlpha',obj.Transparency,'EdgeAlpha',obj.Transparency);
                [~]=patch(ax,'Vertices',ScaledPoints,'Faces',...
                    cell2mat(geometry.BoundaryEdges),'EdgeColor',...
                    obj.EdgeColor,'FaceColor','None');
            end
            axis(ax,'equal')
            %Modified for ticks on shapes that are not on XY plane
            em.MeshGeometry.decoratefigureandaxes(ax, ScaledPoints(:,1), ScaledPoints(:,2), ScaledPoints(:,3), unit, [0 90]);
            box(ax,'on'); grid(ax,'on');  
        end

        % Transformation Methods
        function varargout = translate(obj,vec,varargin)
            %TRANSLATE   Move the shape to a new position
            %
            % OBJ = translate(obj,vec) moves the shape to a new location by
            % using the translation vector.
            %
            % Example
            %
            % % Example1: Create a rectangular shape and translate it to (1,1).
            %
            % r  = shape.Rectangle;
            % r  = translate(r,[1,1,0]);
            %
            % See also <a href="matlab:help shape.Shape2D.rotate">rotate</a>, <a href="matlab:help shape.Shape2D.show">show</a>

            if isempty(obj)
                error(message("antenna:shapeerrors:EmptyShape"));
            end
            validateattributes(vec,{'numeric'},...
                {'nonempty','nonnan','finite','real',...
                'numel',3},'translate','Translation Vector',2);
            if ~isrow(vec)
                vec = vec';
            end
            if nargin==4
                parserObj=inputParser;
                addParameter(parserObj,'Parent',[]);
                parse(parserObj,varargin{:});
                ParentObj=parserObj.Results.Parent;
                checkParentInput(obj,ParentObj);

            else
                ParentObj=[];
            end

            if getHasStructureChanged(obj)
                createGeometry(obj);
            end

            if isequal(vec(3),0)
                % Translate Polygon
                obj.InternalPolyShape = translate(obj.InternalPolyShape,vec(1:2));

                % Update Geometry
                % updateGeometryCache(obj);
                obj.MesherStruct.Geometry.BorderVertices = em.internal.translateshape(obj.MesherStruct.Geometry.BorderVertices',vec)';

                % Center cannot be the property of unsymmetrical shapes the
                % the ReferencePoint will the property to for unsymmetrical PCB
                % shapes to translate w.r.t the reference point.
                if isprop(obj,'ReferencePoint')
                    tempCenter = zeros(1,2);
                    tempCenter(1) = obj.ReferencePoint(1)+vec(1,1) ;
                    tempCenter(2) = obj.ReferencePoint(2)+vec(1,2) ;
                    obj.ReferencePoint = tempCenter;
                else
                    updateCenterProp(obj);
                end
                resetHasStructureChanged(obj);
                % 3. Update Mesh if it exists
                isMeshed = ~isempty(obj.MesherStruct.Mesh.p);

                if isMeshed %--> allow translation of mesh, show
                    obj.MesherStruct.Mesh.p = em.internal.translateshape(obj.MesherStruct.Mesh.p,vec);
                    obj.Points = obj.MesherStruct.Mesh.p;
                end
            else
                translateOutOfXYPlane(obj,vec);
            end
            % Show final structure
            if nargout==0
                 if isempty(ParentObj)
                    show(obj);
                else
                    show(obj,"Parent",ParentObj);
                end
            else
                varargout{1} = obj;
            end

        end

        function varargout = rotate(obj,ang,axispt1,axispt2,varargin)
            %ROTATE   Rotate the shape around an axis by a specific angle
            %
            % OBJ = rotate(obj,ang,axispt1,axispt2) rotates the shape about the
            % specified axis vector defined with the two end points and by the
            % specified angle.
            %
            % Example
            %
            % % Example1: Create a rectangular shape and rotate it by 45
            % % degrees about the z-axis.
            %
            % r  = shape.Rectangle;
            % rotate(r,45,[0,0,0],[0,0,1]);
            %
            % See also <a href="matlab:help shape.Shape2D.translate">translate</a>, <a href="matlab:help shape.Shape2D.scale">scale</a>


            if isempty(obj)
                error(message("antenna:shapeerrors:EmptyShape"));
            end

            nameOfFunction = 'rotate';
            validateattributes(axispt1,{'double'},{'nonempty','nonnan','finite','real',      ...
                                       'numel',3},nameOfFunction,'Axis point1',3);
            validateattributes(axispt2,{'double'},{'nonempty','nonnan','finite','real',      ...
                                       'numel',3},nameOfFunction,'Axis point2',4);
            
            if nargin==6
                parserObj=inputParser;
                addParameter(parserObj,'Parent',[]);
                parse(parserObj,varargin{:});
                ParentObj=parserObj.Results.Parent;
                checkParentInput(obj,ParentObj);
            else
                ParentObj=[];
            end


            if getHasStructureChanged(obj)
                createGeometry(obj);
            end

            tempAng = ang;

            % Rotate
            rvec = axispt2-axispt1;
            if abs(rvec(1))<=sqrt(eps) && abs(rvec(2)) <=sqrt(eps)
                refpoint = axispt1(1:2);
                %Reversing Angle for Reversed Axes
                if rvec(3)<0
                    xyang = -1*ang;
                else
                    xyang = ang;
                end

                obj.InternalPolyShape = rotate(obj.InternalPolyShape,xyang,refpoint);
                obj.MesherStruct.Geometry.BorderVertices = em.internal.rotateshape(obj.MesherStruct.Geometry.BorderVertices',axispt1,axispt2,ang)';

                % 3. Update Mesh if it exists
                isMeshed = ~isempty(obj.MesherStruct.Mesh.p);

                if isMeshed  %--> allow rotation of mesh, show
                    obj.MesherStruct.Mesh.p = em.internal.rotateshape(obj.MesherStruct.Mesh.p,axispt1,axispt2,ang);
                    obj.Points = obj.MesherStruct.Mesh.p;
                    obj.Triangles = obj.MesherStruct.Mesh.t;
                end
                obj.Tilt = [obj.Tilt tempAng];
                tmpAxis = [axispt1;axispt2];
                obj.TiltAxis = [obj.TiltAxis; tmpAxis];
                resetHasStructureChanged(obj);
            else
                rotateAlongXYAxis(obj,ang,axispt1,axispt2);
            end

            % Show final structure
            if nargout==0
                 if isempty(ParentObj)
                    show(obj);
                else
                    show(obj,"Parent",ParentObj);
                end
            else
                varargout{1} = obj;
            end
        end

        function varargout = rotateX(obj,ang,varargin)
            %ROTATEX   Rotate the shape around the x-axis by a specific angle
            %
            % OBJ = rotateX(obj,ang) rotates the shape about the
            % specified axis vector defined with the two end points and by the
            % specified angle.
            %
            % Example
            %
            % % Example1: Create a rectangular shape and rotate it by 45
            % % degrees about the x-axis.
            %
            % r  = shape.Rectangle;
            % rotateX(r,45);
            %
            % See also <a href="matlab:help shape.Shape2D.rotate">rotate</a>, <a href="matlab:help shape.Shape2D.translate">translate</a>


            if isempty(obj)
                error(message("antenna:shapeerrors:EmptyShape"));
            end

            if nargin==4
                parserObj=inputParser;
                addParameter(parserObj,'Parent',[]);
                parse(parserObj,varargin{:});
                ParentObj=parserObj.Results.Parent;
                checkParentInput(obj,ParentObj);
            else
                ParentObj=[];
            end


            if getHasStructureChanged(obj)
                createGeometry(obj);
            end
            rotateAlongXYAxis(obj,ang,[0, 0, 0],[1, 0, 0]);

             % Show final structure
            if nargout==0
                if isempty(ParentObj)
                    show(obj);
                else
                    show(obj,"Parent",ParentObj);
                end
            else
                varargout{1} = obj;
            end

            resetState = getHasStructureChanged(obj);
            if isprop(obj,'Center')
                tmp = em.internal.rotateshape(([obj.Center 0])',[0 0 0],[1 0 0],ang);
                obj.Center = [tmp(1) tmp(2)];
            elseif isprop(obj,'ReferencePoint')
                tmp = em.internal.rotateshape(([obj.ReferencePoint 0])',[1 0 0],[0 0 0],ang);
                obj.ReferencePoint = [tmp(1) tmp(2)];
            end
            obj.MesherStruct.HasStructureChanged = resetState;

        end

        function varargout = rotateY(obj,ang,varargin)
            %ROTATEY   Rotate the shape around the y-axis by a specific angle
            %
            % OBJ = rotateY(obj,ang) rotates the shape about the
            % specified axis vector defined with the two end points and by the
            % specified angle.
            %
            % Example
            %
            % % Example1: Create a rectangular shape and rotate it by 45
            % % degrees about the y-axis.
            %
            % r  = shape.Rectangle;
            % rotateY(r,45);
            %
            % See also <a href="matlab:help shape.Shape2D.rotateX">rotateX</a>, <a href="matlab:help shape.Shape2D.translate">translate</a>


            if isempty(obj)
                error(message("antenna:shapeerrors:EmptyShape"));
            end

            if nargin==4
                parserObj=inputParser;
                addParameter(parserObj,'Parent',[]);
                parse(parserObj,varargin{:});
                ParentObj=parserObj.Results.Parent;
                checkParentInput(obj,ParentObj);
            else
                ParentObj=[];
            end


            if getHasStructureChanged(obj)
                createGeometry(obj);
            end
            rotateAlongXYAxis(obj,ang,[0 0 0],[0 1 0]);
            % Show final structure
            if nargout==0
                if isempty(ParentObj)
                    show(obj);
                else
                    show(obj,"Parent",ParentObj);
                end
            else
                varargout{1} = obj;
            end

            resetState = getHasStructureChanged(obj);
            if isprop(obj,'Center')
                tmp = em.internal.rotateshape(([obj.Center 0])',[0 0 0],[0 1 0],ang);
                obj.Center = [tmp(1) tmp(2)];
            elseif isprop(obj,'ReferencePoint')
                tmp = em.internal.rotateshape(([obj.ReferencePoint 0])',[0 0 0],[0 1 0],ang);
                obj.ReferencePoint = [tmp(1) tmp(2)];
            end
            obj.MesherStruct.HasStructureChanged = resetState;

        end

        function varargout = rotateZ(obj,ang,varargin)
            %ROTATEZ   Rotate the shape around the z-axis by a specific angle
            %
            % OBJ = rotateZ(obj,ang) rotates the shape about the
            % specified axis vector defined with the two end points and by the
            % specified angle.
            %
            % Example
            %
            % % Example1: Create a rectangular shape and rotate it by 45
            % % degrees about the z-axis.
            %
            % r  = shape.Rectangle;
            % rotateZ(r,45);
            %
            % See also <a href="matlab:help shape.Shape2D.rotateX">rotateX</a>, <a href="matlab:help shape.Shape2D.translate">translate</a>


            if isempty(obj)
                error(message("antenna:shapeerrors:EmptyShape"));
            end
            if nargin==4
                parserObj=inputParser;
                addParameter(parserObj,'Parent',[]);
                parse(parserObj,varargin{:});
                ParentObj=parserObj.Results.Parent;
                checkParentInput(obj,ParentObj);
            else
                ParentObj=[];
            end

            % The ReferencePoint will be the property for unsymmetrical
            % shapes to rotate the shape in reference to the points
            obj = rotate(obj,ang,[0, 0, 0],[0, 0, 1]);

             % Show final structure
            if nargout==0
                if isempty(ParentObj)
                    show(obj);
                else
                    show(obj,"Parent",ParentObj);
                end
            else
                varargout{1} = obj;
            end

        end

        function varargout = scale(obj,scaling,varargin)
            %SCALE   Change the size of the shape by a fixed amount
            %
            % OBJ = scale(obj,scaling) scales the shape by a constant factor
            %
            % Example
            %
            % % Example1: Create a rectangular shape and shrink it by 50%.
            %
            % r  = shape.Rectangle;
            % scale(r,0.5);
            %
            % See also <a href="matlab:help shape.Shape2D.rotateZ">rotateZ</a>, <a href="matlab:help shape.Shape2D.translate">translate</a>

            if isempty(obj)
                error(message("antenna:shapeerrors:EmptyShape"));
            end
            if nargin == 2
                RefPoint = [];
            end
            if nargin==3||nargin==5
                RefPoint = varargin{1};
                 validateattributes(RefPoint,{'numeric'},...
                    {'nonempty','nonnan','finite','real',...
                    'numel',2},'scale','RefPoint',3);
            end
            if nargin==4 
                RefPoint=[];
                parserObj=inputParser;
                addParameter(parserObj,'Parent',[]);
                parse(parserObj,varargin{1:2});
                ParentObj=parserObj.Results.Parent;
            elseif nargin ==5
                parserObj=inputParser;
                addParameter(parserObj,'Parent',[]);
                parse(parserObj,varargin{2:3});
                ParentObj=parserObj.Results.Parent;
            else
                ParentObj=[];
            end
           
            checkParentInput(obj,ParentObj);
            
            if getHasStructureChanged(obj)
                createGeometry(obj);
            end
            if isa(obj,'shape.Circle')&& ~(isscalar(scaling))
                error(message('antenna:shapeerrors:OnlyScalar', ...
                    'scaling for shape.Circle'));
            end
            % Scale
            % centroid cannot hold good on unsymmetrical shapes. The
            % the ReferencePoint will the property to for unsymmetrical RF PCB
            % shapes shape the shapes based on the reference point of the shapes.
            S = scaling;
            if isprop(obj,'ReferencePoint')
                RefPoint_x = obj.ReferencePoint(1);
                RefPoint_y = obj.ReferencePoint(2);
                obj.InternalPolyShape = scale(obj.InternalPolyShape,S,[RefPoint_x,RefPoint_y]);
            else
                [xcenter ,ycenter]  = centroid(obj.InternalPolyShape);
                if ~isempty(RefPoint)
                    % Check if the shape is symmetrical about x and y axis.
                    % The scale operation, with respect to the centroid is valid for
                    % all the symmetrical shapes as the centroid of the
                    % shape is the center. Hence the reference point has to
                    % be ignored if the shape is symmetrical and consider
                    % if it is not.
                    V = getShapeVertices(obj);
                    V = unique(V,'rows','stable');
                    A = V(:,1:2);
                    % sort the y points
                    A2y = unique((sort(A(:,2),'ascend')));
                    % Find x points for the sorted y points
                    x1 = cell(1,size(A2y,1));
                    x = [];
                    for i = 1:size(A2y,1)
                        x1{i} = find(abs(A(:,2)-A2y(i)) <= eps(max(A(:,2),A2y(i))));
                        %x1{i} = find(A(:,2)==A2y(i));
                        x(end+1:end+size(x1{i},1),:)     = x1{i};
                    end
                    x2 = unique(x,'row','stable');
                    Ax = A(x2,1);
                    Ay = sort(A(:,2),'ascend');
                    % Vertices after sorting
                    A2 = [Ax,Ay];
                    % check the distances from the centroid to all the points
                    distances = round(sqrt((A2(:,1)-xcenter).^2 + (A2(:,2)-ycenter).^2),12);
                    n = numel(A2)/2;
                    if mod(n,2)==1
                        n1 = round(floor(n/2),12);
                    else 
                        n1 = n/2;
                    end
                    % Distance from the top to the mid points and from the bottom to the mid
                    % points should be same if the structure is symmetry 
                    d1 = distances(1:n1);
                    d1 = sort(d1,'ascend');
                    d2 = flip(distances);
                    d3 = d2(1:n1);
                    d3 = sort(d3,'ascend');
                    if all(d1==d3)||isprop(obj,'Center')
                        warning(message('antenna:shapeerrors:IgnoreRefPoint'))
                        center_x =  xcenter;
                        center_y =  ycenter;
                    else
                        center_x =  RefPoint(1);
                        center_y =  RefPoint(2);
                    end
                else
                    center_x =  xcenter;
                    center_y =  ycenter;
                end
                obj.InternalPolyShape = scale(obj.InternalPolyShape,S,[center_x center_y]);
            end
            
            updateShapeDimension(obj,S);
            resetHasStructureChanged(obj);
            % 3. Update Mesh if it exists
            isMeshed = ~isempty(obj.MesherStruct.Mesh.p);

            if isMeshed
                obj.MesherStruct.Mesh.p = (S*obj.MesherStruct.Mesh.p);
                obj.Points = obj.MesherStruct.Mesh.p;
            end

            if nargout==0
                if isempty(ParentObj)
                    show(obj);
                else
                    show(obj,"Parent",ParentObj);
                end
            else
                varargout{1} = obj;
            end
        end

        % Calculation Methods

        function A = area(obj)
            %AREA   Find the area of the shape in sq.m
            %
            % A = area(obj) computes the area of the shape
            %
            % Example
            %
            % % Example1: Create a rectangular shape and find its area.
            %
            % r  = shape.Rectangle;
            % A = area(r);
            %
            % See also <a href="matlab:help shape.Shape2D.mesh">mesh</a>, <a href="matlab:help shape.Shape2D.translate">translate</a>

            if isempty(obj)
                error(message("antenna:shapeerrors:EmptyShape"));
            end

            if getHasStructureChanged(obj)
                createGeometry(obj);
            end
            A = area(obj.InternalPolyShape);
        end

        function [x,y] = centroid(obj)
            %CENTROID   Find the centroid of the shape in sq.m
            %
            % [x,y] = centroid(obj) computes the area of the shape
            %
            % Example
            %
            % % Example1: Create a rectangular shape and find its area.
            %
            % r  = shape.Rectangle;
            % [x,y] = centroid(r);
            %
            % See also <a href="matlab:help shape.Shape2D.area">area</a>, <a href="matlab:help shape.Shape2D.translate">translate</a>

            if getHasStructureChanged(obj)
                createGeometry(obj);
            end
            [x,y] = centroid(obj.InternalPolyShape);
        end

        function [in,on] = contains(obj,x,y)
            if getHasStructureChanged(obj)
                createGeometry(obj);
            end
            [in,on] = isinterior(obj.InternalPolyShape,x,y);
            %             [in,~] = isinterior(obj.InternalPolyShape,x,y);
        end

        % Clean up methods

        function s = removeSlivers(obj,slivertol)
            %   S = REMOVESLIVERS(OBJ,SLIVERTOL) removes sliver outliers
            %   from boundary of shape
            %
            % Example
            %
            % % Example1: Remove the sliver-like protrusion from a
            % rectangular shape
            %
            % a = shape.Rectangle;
            % b = shape.Rectangle;
            % b.Length = 1e-7;
            % b.Center = [0,0.75];
            % c = a+b;
            % d = removeSlivers(c,1e-6);
            %
            % See also <a href="matlab:help shape.Shape2D.plus">plus</a>, <a href="matlab:help shape.Shape2D.removeHoles">area</a>


            layerpoly = obj.InternalPolyShape;
            fixedpoly = rmslivers(layerpoly,slivertol);
            fixedobject = shape.Polygon('Vertices',fixedpoly.Vertices);
            s = fixedobject;
        end

        function s = removeHoles(obj,holetol)
            % S = REMOVEHOLES(OBJ,HOLETOL) removes holes with area less
            % than the provided tolerance from the shape
            %
            % Example
            %
            % % Example1: Remove the hole from a rectangular shape of area
            % less than 1e-6
            %
            % a = shape.Rectangle;
            % b = shape.Rectangle;
            % b.Length = 1e-7;
            % b.Width = 0.5
            % c = a-b;
            % d = removeHoles(c,1e-6);
            %
            % See also <a href="matlab:help shape.Shape2D.minus">minus</a>, <a href="matlab:help shape.Shape2D.removeSlivers">removeSlivers</a>

            layerpoly = obj.InternalPolyShape;
            if layerpoly.NumHoles>0
                h = holes(layerpoly);
                % Test each hole against the hole tol area
                holeareas = area(h);
                id = find(holeareas<holetol);
                % Eliminate vertices belonging to holes smaller than hole tol
                fixedpoly = rmboundary(layerpoly,id+1); % Outermost boundary is the first..
                fixedobject = shape.Polygon('Vertices',fixedpoly.Vertices);
                s = fixedobject;
            end
        end

        % Access Shape Vertices

        function v = getShapeVertices(obj)
            if getHasStructureChanged(obj)
                createGeometry(obj);
            end
            v = obj.ShapeVertices;
        end

    end

    methods(Hidden)

        function [edgeLength,growthRate] = calculateMeshParams(obj,varargin)
            A       = area(obj);
            lambda  = varargin{1};
            N       = A/(lambda^2);
            frac    = 10; %varargin{2};
            edgeLength       = sqrt((4*A)/(frac*N*100*sqrt(3)));
            growthRate = 1.7;
        end

        function meshGenerator(obj,varargin)
            % Check if structure has changed and regenerate mesh if needed
            %
            if getHasStructureChanged(obj)
                createGeometry(obj);
            end

            % Check if Shape is still on XY
            if ~isShapeOnXY(obj)
                error(message('antenna:shapeerrors:Unsupported', ...
                    'Boolean operation','shapes not on XY-plane at z-coordinate of 0. Check the shape'));
            end

            % Mesh structure
            if numel(varargin)>1
                edgeLength  = varargin{1};
                growthRate  = varargin{2};
                setMeshEdgeLength(obj, edgeLength);
                setMeshGrowthRate(obj, growthRate);
                if numel(varargin) >2
                    minEdgeLength = varargin{3};
                    targetEdgeLength = varargin{4};
                    minFeatureEdgeLength = varargin{5};
                    setMeshMinContourEdgeLength(obj,minEdgeLength);
                    setTargetInteriorEdgeLength(obj,targetEdgeLength);
                    setSubDomainEdgeLength(obj,minFeatureEdgeLength);
                end
            else
                edgeLength  = getMeshEdgeLength(obj);
                growthRate  = getMeshGrowthRate(obj);
                minEdgeLength = getMinContourEdgeLength(obj);
            end

            [p,t] = getInitialMesh(obj);
%             FeedsAndVias = obj.MesherStruct.Mesh.PartMesh.Others.p;
            setFillAndHolePolygons(obj);
            if ~isempty(obj.FeaturePolygons)
                [p,t] = fillFeatureHoleInMesh(obj,p,t);
            end

            % Check mesh for triangles with small areas
            A = em.internal.meshprinting.meshareas(p, t);
            % Find triangles with areas approaching machine epsilon and remove
            Epsilon = 1e-12; %sqrt(eps);
            id = A<Epsilon;
            t(id,:) = [];
            FeedsAndVias=[];

            [ps,es,ts] = em.MeshGeometry.remesher(p,t,edgeLength,growthRate,minEdgeLength,...
                targetEdgeLength,minFeatureEdgeLength,  ...
                FeedsAndVias,[]);
            ts(4,:) = 0;

            % Check if less than 10 triangles are generated and if so error
            %             trilimit = 5;
            %             if size(ts,2) < trilimit
            %                 error(message('antenna:shapeerrors:InsufficientTrianglesInMesh',num2str(trilimit)));
            %             end
            
            obj.Edges = es;
            T           = [];
            EpsilonR    = [];
            LossTangent = [];
            Mesh        = em.internal.makeMeshStructure(ps,ts,T,EpsilonR,...
                LossTangent);

            saveMesh(obj,Mesh);
        end

        function s = saveobj(obj)           %---SAVEOBJ IMPLEMENTATION----%
            obj.MesherStruct.CacheFlag = obj.MesherStruct.HasStructureChanged;
            s = obj;
        end

        function setFeaturePolygons(obj,h)
            obj.FeaturePolygons = h;
        end

        function setFillAndHolePolygons(obj)
            polyObj     = obj.InternalPolyShape;
            % Add the feature polygons back in
            if ~isempty(obj.FeaturePolygons)
                for i = 1:numel(obj.FeaturePolygons)
                    polyObj = union(polyObj,obj.FeaturePolygons{i}.InternalPolyShape);
                end
            end
            fillPolys   = regions(rmholes(polyObj));
            obj.FillPolygons = arrayfun(@(x) x.Vertices,fillPolys,'UniformOutput',false)';
            holePolys   = holes(polyObj);
            % If there are Non-zero hole polygons
            if ~isempty(holePolys)
                obj.HolePolygons = arrayfun(@(x) x.Vertices,holePolys,'UniformOutput',false)';
            end

            %             numFeatureHoles  = numel(featurePolys);
            %             if ~isequal(numFeatureHoles,0)&&~isempty(holePolys)
            %                 featureHoleId = nan.*ones(1,numFeatureHoles);
            %                 for i = 1:numFeatureHoles
            %                     tempShape = arrayfun(@(x) (subtract(x,featurePolys{i}.InternalPolyShape)),holePolys);
            %                     isFeatureHole = arrayfun(@(y) isempty(y.Vertices),tempShape);
            %                     featureHoleId(i) = find(isFeatureHole);
            %                 end
            %
            %                 % If there are feature polygons
            %                 if ~isempty(featureHoleId)
            %                     % Remove from list of hole polygons
            %                     holePolys(featureHoleId) = [];
            %                 end
            %             end

        end

        function updateShape(obj,propVal)
            obj.InternalPolyShape.Vertices = propVal(:,1:2);
        end

       
        function tf = hasFeatureVerticesOnLayer(obj,f)
            layerPoly   = obj.InternalPolyShape;
            featurePoly = f.InternalPolyShape;
            %fv = intersect(layerPoly.Vertices,featurePoly.Vertices,'rows','stable');
            lia = ismembertol(layerPoly.Vertices,featurePoly.Vertices,'ByRows',true);
            fv = layerPoly.Vertices(lia,:);
            tf = true;
            if ~isequal(size(fv,1),size(featurePoly.Vertices,1))
                tf = false;
            end
        end

        function tf = hasGeometry(obj)
            geometry        = getGeometry(obj);
            if isempty(geometry.BorderVertices)
                tf = false;
            else
                tf = true;
            end
        end
        %
        function [p,t] = getInitialMesh(obj)
            TR = triangulation(obj.InternalPolyShape);
            p = TR.Points;
            p(:,3) = 0;
            t = TR.ConnectivityList;
            t(:,4) = 0;
        end

        function [p,t] = fillFeatureHoleInMesh(obj,p,t)
            % Find number of feature holes to add back into mesh
            numFeatureHoles = numel(obj.FeaturePolygons);
            % Get layer triangulation
            TRPoly = triangulation(obj.InternalPolyShape);
            holeTriangles = [];
            for i =1:numFeatureHoles
                trhole = triangulation(obj.FeaturePolygons{i}.InternalPolyShape);
                % Find hole points in layer
                ia = find(ismembertol(TRPoly.Points,trhole.Points,sqrt(eps),'DataScale',1,'ByRows',true));
                % Use as index map into connectivity for building triangles
                holeTriangles = [holeTriangles;ia(trhole.ConnectivityList)]; %#ok<AGROW>
            end
            holeTriangles(:,4) = 0;
            % Repair triangulation
            t = [t;holeTriangles];
        end

    end

    methods(Access={?shape.CADEngine})
         function tf = hasNoSlivers(obj)
            tf = true;
            if obj.InternalPolyShape.NumRegions>1
                r = regions(obj.InternalPolyShape);
                A = area(r);
                if any(A<obj.SliverTol)
                    tf = false;
                end
            end
         end

         function tf = isShapeEmpty(obj)
             tf = isempty(obj.InternalPolyShape.Vertices);
         end

         function   assignMultiShapeIDFor2D(obj,s1,s2)
           % Set the multishapeid for first shape,if it is empty set it as 1,else
            % increase the number
            tr1=triangulation(s1.InternalPolyShape);
            tr2 = triangulation(s2.InternalPolyShape);
            ints = intersect(s1.InternalPolyShape,s2.InternalPolyShape);
            if isempty(ints.Vertices)
                if isempty(s1.MultiShapeID)
                    polyg1=ones(size(tr1.ConnectivityList,1),1);
                    s1.MultiShapeID=polyg1;
                else
                    polyg1=s1.MultiShapeID;
                end

            % Set the multishapeid for shape2,if it is empty,increase the number from
            % max of shape1
            if isempty(s2.MultiShapeID)
                if isempty(s1.MultiShapeID)
                    polyg2=ones(size(tr2.ConnectivityList,1),1)*2;
                    s2.MultiShapeID=polyg2;
                else
                    polyg2num=max(s1.MultiShapeID)+1;
                    s2.MultiShapeID=ones(size(tr2.ConnectivityList,1),1)*polyg2num;
                    polyg2=s2.MultiShapeID;
                end
            else
                %     trnew1(:,4)=s2.MultiShapeID;
                polyg2=max(s1.MultiShapeID)+s2.MultiShapeID;
                s2.MultiShapeID=polyg2;
            end
            obj.MultiShapeID=[polyg1;polyg2];
            assignColorFor2D(obj,s1,s2,1);
            addcolorproperties(obj);

            else
%              obj.MultiShapeID = ones(size(cell2mat(obj.Polygons),1),1);
             obj.MultiShapeID=[];
             assignColorFor2D(obj,s1,s2,0);


            end
         end
         function    assignColorFor2D(obj,s1,s2,fl)
            [~,~]=getshapeparameters(obj,s1,s2);
            % Assign color id
            obj.ColorID={};
            obj.ColorID(end+1:end+size(s1.ColorID,1),:)=s1.ColorID;
            if fl
                if size(s2.ColorID,1)==1
                    siz=size(s1.ColorID,1);
                    obj.ColorID(siz+1:siz+size(s2.ColorID,1),:)=s2.ColorID ;
                end
            end
            s1.MultiShapeID=[];s2.MultiShapeID=[];
        end

    end

    methods(Access=protected)

        function rotvert = rotateAngVert(~,vert, angle, NumSegments)
            % Generate the total vertices based on step angle .
            stepangle = angle/NumSegments;
            rotvert = zeros(NumSegments*size(vert,1),3);
            n = size(vert,1);
            for i = 0:NumSegments
                rotvert(i*n +1 : i*n + n,:) = transpose(em.internal.rotateshape(vert',[0,-1,0],[0,1,0],stepangle*i));
            end
        end

        function transVert = translateVert(shapin,vert, height, NumSegments, varargin)
            % This function translates the vertices of 2d shape step by step based on
            % the height and numsegments
            twist = 0;
            axis = mean(vert);
            twistradius = 0;
            if ~isempty(varargin)
                twist = varargin{1};
                twistradius =varargin{2};
                Scale=varargin{3};
            end
            % Taking the effect of twist radius
            axis(1) = twistradius + axis(1);
            %points changing with scalar and vector inputs
            if isscalar(Scale)
                twradnew1=linspace(1,Scale,NumSegments+1);
                twradnew2=linspace(1,Scale,NumSegments+1);

            else
                twradnew1=linspace(1,Scale(1),NumSegments+1);
                twradnew2=linspace(1,Scale(2),NumSegments+1);
            end
            stepheight = height/NumSegments;
            transVert = zeros(NumSegments*size(vert,1),3);
            n = size(vert,1);
            %      if ~isscalar(Scale)
            %          vert=vert*Scale(1);
            %      end
            vertnew=[];
            for i = 0:NumSegments
                % Translate the vertices based on upon the twist
                if twist*i == 0
                    if ~(Scale(1)==1)
                        vertnew(:,1)=vert(:,1)*twradnew1(i+1);
                        vertnew(:,2)=vert(:,2)*twradnew2(i+1);
                        vertnew(:,3)=vert(:,3);
                        transVert(i*n +1 : i*n + n,:) = rotatevert(shapin,vertnew,twist*i,axis) +[0,0, stepheight*i];
                    else
                        transVert(i*n +1 : i*n + n,:) = vert +[0,0, stepheight*i];
                    end
                else
                    vertnew(:,1)=vert(:,1)*twradnew1(i+1);
                    vertnew(:,2)=vert(:,2)*twradnew2(i+1);
                    vertnew(:,3)=vert(:,3) ;
                    transVert(i*n +1 : i*n + n,:) = rotatevert(shapin,vertnew,twist*i,axis) +[0,0, stepheight*i];
                end
            end
        end

        function vert = rotatevert(~,inp,angle,axis)
            vert =  transpose(em.internal.rotateshape(inp',[axis(1:2),-1],[axis(1:2),1],angle));
        end
        function conn = genTriForExtrude(~,numvert,NumSegments, rotateflag)
            % Generate the connections
            vertindices = 1:numvert;
            vertindices = vertindices(:);
            conn = [];
            if rotateflag
                firstindex = vertindices;
                secondindex = vertindices+ numvert*(NumSegments-1);%numseg for linear
                conn = [conn;...
                    connectIndicesTriangulation(firstindex,secondindex)];
            end
            % As numsegments increase,connections increases
            for i = 1:NumSegments
                firstindex = vertindices + numvert*(i-1);
                secondindex = vertindices+ numvert*(i);
                %     conn = [conn;...
                %         connectIndicesTriangulation(firstindex,secondindex)];
                cotri=[firstindex, secondindex,circshift(firstindex,-1);...
                    circshift(firstindex,-1),secondindex,circshift(secondindex,-1)];
                conn(end+1:end+size(cotri,1),:)=cotri;
            end
        end

          


 
        function rotateAlongXYAxis(obj,ang,axispt1,axispt2)

            % 1. Update BorderVertices
            rotateVertices(obj,ang,axispt1,axispt2);

            % 2. Update internal object
            if ~any(obj.MesherStruct.Geometry.BorderVertices(:,3)>0)&&  ...
                    ~any(obj.MesherStruct.Geometry.BorderVertices(:,3)<0)
                 nanid = find(isnan(obj.InternalPolyShape.Vertices(:,1)));
                v = obj.InternalPolyShape.Vertices;
                v = [v, zeros(size(v,1),1)];
                if ~isempty(nanid)
                    % Dummy fill
                    v(nanid,:) = repmat([0 0 0],numel(nanid),1);
                    % rotate
                    vr = em.internal.rotateshape(v',axispt1,axispt2,ang)';
                    % restore nan
                    vr(nanid,:) = nan;
                else
                    vr = em.internal.rotateshape(v',axispt1,axispt2,ang)';
                end

                obj.InternalPolyShape = polyshape(vr(:,1:2),'Simplify',false,'KeepCollinearPoints',true);
                obj.Tilt = [obj.Tilt ang];
                tmpAxis = [axispt1;axispt2];
                obj.TiltAxis = [obj.TiltAxis; tmpAxis];
            end

            % 2. Update Mesh if it exists
            isMeshed = ~isempty(obj.MesherStruct.Mesh.p);

            if isMeshed  %--> allow rotation of mesh, show
                obj.MesherStruct.Mesh.p = em.internal.rotateshape(obj.MesherStruct.Mesh.p,axispt1,axispt2,ang);
                obj.Points = obj.MesherStruct.Mesh.p;
            end

        end

        function rotateVertices(obj,ang,axispt1,axispt2)
            Geometry                = getGeometry(obj);
            obj.MesherStruct.Geometry.BorderVertices = em.internal.rotateshape(Geometry.BorderVertices'...
                ,axispt1,axispt2,ang)';
        end

        function translateOutOfXYPlane(obj,vec)
            % 1. Update BorderVertices
            translateVertices(obj,vec)

            % 2. Update Mesh if it exists
            isMeshed = ~isempty(obj.MesherStruct.Mesh.p);

            if isMeshed  %--> allow rotation of mesh, show
                obj.MesherStruct.Mesh.p = em.internal.translateshape(obj.MesherStruct.Mesh.p,vec);
                obj.Points = obj.MesherStruct.Mesh.p;
            end

        end
        % %
        function translateVertices(obj,vec)
            Geometry                = getGeometry(obj);
            obj.MesherStruct.Geometry.BorderVertices = em.internal.translateshape(Geometry.BorderVertices'...
                ,vec)';
        end

 
        function updateCenterProp(obj)
            tmpCenter = zeros(1,2);
            if isprop(obj,'Center')
                [tmpCenter(1),tmpCenter(2)] = centroid(obj);
                idxLTeps = abs(tmpCenter) < eps;
                tmpCenter(idxLTeps) = 0;
                obj.Center = tmpCenter;
            end
        end

        function updateShapeDimension(obj, ~)
            if isa(obj,'shape.Polygon')
                return;
            end
        end
    end

    methods(Hidden)
         function tf = isShapeOnXY(obj)
            % Check if coordinates are not on x-y plane  and returns the
            % appropriate flag
            tf = true;
            polyBorderVerts = obj.getGeometry.BorderVertices;
            if ~isempty(polyBorderVerts)
                isNanIds        = isnan(polyBorderVerts(:,1));
                nanIds          = isNanIds==1;
                % Remove nans temporarily
                polyBorderVerts(nanIds,:) = [];
                if any(abs(polyBorderVerts(:,3))>=sqrt(eps))
                    tf = false;
                end
            end
        end


    end

    methods
        function m = meshconfig(obj,mode)
            % MESHCONFIG Change the mesh mode for shape
            %
            % m = meshconfig(obj,mode), changes the meshing mode for the shape
            % according to the string input mode. The two modes provided are
            % manual and auto. In auto, the highest frequency provided in any analysis
            % is automatically used for meshing the structure. The manual mode provides
            % more control on the mesh for the user. This is done by specifying the
            % maximum edge length in the <a href="matlab:help
            % shape.Shape2D.mesh">mesh</a>.
            %
            % % Examples:
            %
            % % Example 1: Create a rectangle and choose to manually mesh it
            % with a maximum edge length of .015.
            %
            % m = shape.Rectangle;
            % meshconfig(m,'manual')
            % f = mesh(m,MaxEdgeLength = 0.15)
            %
            % See also <a href="matlab:help shape.Shape2D.mesh">mesh</a>,<a href="matlab:help shape.Shape2D.show">show</a>


            m = protectedmeshconfig(obj,mode);
        end

    end

    methods(Hidden)
        function rObj = superLoadShape(obj,s)
                obj.FillPolygons = s.FillPolygons;
                obj.HolePolygons = s.HolePolygons;
                if isfield(s,'FeaturePolygons') || isprop(s,'FeaturePolygons')
                    obj.FeaturePolygons = s.FeaturePolygons;
                end
                if isfield(s,'InternalPolyShape')||isprop(s,'InternalPolyShape')
                    obj.InternalPolyShape = s.InternalPolyShape;
                end
                obj.Points = s.Points;
                obj.Triangles =  s.Triangles;
                obj.Edges = s.Edges;
                obj.MultiShapeID =s.MultiShapeID;
                obj.Metal=s.Metal;
                obj.Color = s.Color;
                obj.Transparency =s.Transparency;
                obj.EdgeColor = s.EdgeColor;
                obj.MaterialID = s.MaterialID;
                obj.TetID   = s.TetID;
                obj.ColorID = s.ColorID;
                obj.HasMaterialID = s.HasMaterialID;
                obj.HasTetID   = s.HasTetID;
                rObj = obj;
           
        end
    end

    methods(Static,Hidden)

        function B = findBoundary(TR)
            % [~,vid] = boundaryshape(TR);
            vid = freeBoundary(TR);
            nanid = find(isnan(vid));
            if ~isempty(nanid)
                startid = 1;
                bo = cell(1,numel(nanid));
                for i = 1:numel(nanid)
                    stopid = nanid(i)-1;
                    bo{i} = [vid(startid:stopid);vid(startid)];
                    startid = nanid(i) + 1;
                end

                bo{i+1} = [vid(startid:end);vid(startid)];
                N = cellfun(@(x) (numel(x)),bo);
                Nmax = max(N);
                B = nan.*ones(Nmax,numel(N));
                for i = 1:numel(N)
                    B(1:N(i),i) = bo{i};
                end
            else
                B = vid;
            end
            B = B;
        end
        function r = loadobj(obj)
            % checkFeatureFlag(obj);
            r = obj;
            r.MesherStruct.HasStructureChanged = obj.MesherStruct.CacheFlag;
        end
    end
end
