diff --git a/.github/workflows/nightly-doc-build.yml b/.github/workflows/nightly-doc-build.yml
index 714468dc97..870d41580c 100644
--- a/.github/workflows/nightly-doc-build.yml
+++ b/.github/workflows/nightly-doc-build.yml
@@ -56,7 +56,6 @@ jobs:
 
       - name: Deploy
         uses: JamesIves/github-pages-deploy-action@3.7.1
-        if: startsWith(github.ref, 'refs/tags/')
         with:
           repository-name: pyansys/pymapdl-dev-docs
           token: ${{ secrets.PYANSYS_CI_BOT_TOKEN }}
diff --git a/ansys/mapdl/core/_commands/preproc/elements.py b/ansys/mapdl/core/_commands/preproc/elements.py
index 9c1b100b52..b54082402d 100644
--- a/ansys/mapdl/core/_commands/preproc/elements.py
+++ b/ansys/mapdl/core/_commands/preproc/elements.py
@@ -1201,7 +1201,8 @@ def eplot(self, show_node_numbering=False, vtk=None, **kwargs):
                 labels = [{'points': esurf.points,
                            'labels': esurf['ansys_node_num']}]
 
-            return general_plotter([{'mesh': esurf}],
+            return general_plotter([{'mesh': esurf,
+                                     'style': kwargs.pop('style', 'surface')}],
                                    [],
                                    labels,
                                    **kwargs)
diff --git a/ansys/mapdl/core/inline_functions/inline_functions.py b/ansys/mapdl/core/inline_functions/inline_functions.py
index 5b08a3543c..cc7e6f5a14 100644
--- a/ansys/mapdl/core/inline_functions/inline_functions.py
+++ b/ansys/mapdl/core/inline_functions/inline_functions.py
@@ -4,6 +4,8 @@
     _NextSelectedEntityQueries
 from .line_queries import _LineFractionCoordinateQueries, \
     _LineFractionSlopeQueries
+from .normals_queries import _NodeNormalQueries, _KeypointNormalQueries
+from .nearest_queries import _EntityNearestEntityQueries
 
 
 class Query(_ComponentQueries,
@@ -12,7 +14,10 @@ class Query(_ComponentQueries,
             _SelectionStatusQueries,
             _LineFractionCoordinateQueries,
             _LineFractionSlopeQueries,
-            _NextSelectedEntityQueries):
+            _NextSelectedEntityQueries,
+            _NodeNormalQueries,
+            _KeypointNormalQueries,
+            _EntityNearestEntityQueries):
     """Class containing all the inline functions of APDL.
 
     Most of the results of these methods are shortcuts for specific
@@ -53,6 +58,9 @@ class Query(_ComponentQueries,
     - ``arnext(a)`` - get the next selected area with a number greater than `a`.
     - ``elnext(e)`` - get the next selected element with a number greater than `e`.
     - ``vlnext(v)`` - get the next selected volume with a number greater than `v`.
+    - ``nnear(n)`` - get the selected node nearest node `n`.
+    - ``knear(k)`` - get the selected keypoint nearest keypoint `k`.
+    - ``enearn(n)`` - get  the selected element nearest node `n`.
     - ``node(x, y, z)`` - get the node closest to coordinate (x, y, z)
     - ``kp(x, y, z)`` - get the keypoint closest to coordinate (x, y, z)
 
diff --git a/ansys/mapdl/core/inline_functions/nearest_queries.py b/ansys/mapdl/core/inline_functions/nearest_queries.py
new file mode 100644
index 0000000000..bb741e431e
--- /dev/null
+++ b/ansys/mapdl/core/inline_functions/nearest_queries.py
@@ -0,0 +1,114 @@
+
+class _EntityNearestEntityQueries:
+    _mapdl = None
+
+    def nnear(self, n: int) -> int:
+        """Returns the selected node nearest node `n`.
+
+        Parameters
+        ----------
+        n : int
+            Node number
+
+        Returns
+        -------
+        int
+            Node number
+
+        Examples
+        --------
+        In this example we construct a solid box and mesh it. Then we
+        use the ``Query`` method ``node`` to find the node closest to
+        the centre of the block and finally use ``nnear`` to find the
+        node nearest to that.
+
+        >>> from ansys.mapdl.core import launch_mapdl
+        >>> mapdl = launch_mapdl()
+        >>> mapdl.prep7()
+        >>> mapdl.et(1, 'SOLID5')
+        >>> mapdl.block(0, 10, 0, 10, 0, 10)
+        >>> mapdl.esize(3)
+        >>> mapdl.vmesh('ALL')
+        >>> q = mapdl.query()
+        >>> node_number = q.node(5., 5., 5.)
+        >>> nearest_node = q.nnear(node_number)
+        >>> node_number, nearest_node
+        (112, 103)
+        """
+        response = self._mapdl.run(f'_=NNEAR({n})')
+        integer = self._parse_parameter_integer_response(response)
+        return integer
+
+    def knear(self, k: int) -> int:
+        """Returns the selected keypoint nearest keypoint `k`.
+
+        Parameters
+        ----------
+        k : int
+            Keypoint number
+
+        Returns
+        -------
+        int
+            Keypoint number
+
+        Examples
+        --------
+        In this example we construct two keypoints and then verify that
+        the nearest keypoint to one is the other.
+
+        >>> from ansys.mapdl.core import launch_mapdl
+        >>> mapdl = launch_mapdl()
+        >>> mapdl.prep7()
+        >>> k1 = mapdl.k(1, 0, 0, 0)
+        >>> k2 = mapdl.k(2, 1, 1, 1)
+        >>> q = mapdl.query()
+        >>> q.knear(k1)
+        1
+
+        >>> q.knear(k1) == k2
+        True
+        """
+        response = self._mapdl.run(f'_=KNEAR({k})')
+        integer = self._parse_parameter_integer_response(response)
+        return integer
+
+    def enearn(self, n: int) -> int:
+        """Returns the selected element nearest node `n`.
+
+        The element position is calculated from the selected nodes.
+
+        Parameters
+        ----------
+        n : int
+            Node number
+
+        Returns
+        -------
+        int
+            Element number
+
+        Examples
+        --------
+        In this example we construct a solid box and mesh it. Then we
+        use the ``Query`` method ``node`` to find the node closest to
+        the centre of the block and finally use ``enearn`` to find the
+        element nearest to that node.
+
+        >>> from ansys.mapdl.core import launch_mapdl
+        >>> mapdl = launch_mapdl()
+        >>> mapdl.prep7()
+        >>> mapdl.et(1, 'SOLID5')
+        >>> mapdl.block(0, 10, 0, 10, 0, 10)
+        >>> mapdl.esize(3)
+        >>> mapdl.vmesh('ALL')
+        >>> q = mapdl.query()
+        >>> node_number = q.node(5., 5., 5.)
+        >>> nearest_element = q.enearn(node_number)
+        >>> node_number, nearest_element
+        (112, 22)
+        """
+        response = self._mapdl.run(f'_=ENEARN({n})')
+        integer = self._parse_parameter_integer_response(response)
+        return integer
+
diff --git a/ansys/mapdl/core/inline_functions/normals_queries.py b/ansys/mapdl/core/inline_functions/normals_queries.py
new file mode 100644
index 0000000000..1a19301ced
--- /dev/null
+++ b/ansys/mapdl/core/inline_functions/normals_queries.py
@@ -0,0 +1,271 @@
+from .core import _ParameterParsing
+
+
+class _NodeNormalQueries(_ParameterParsing):
+    _mapdl = None
+
+    def normnx(self, n1: int, n2: int, n3: int) -> float:
+        """X-direction cosine of the normal to the plane containing the given nodes.
+
+        X-direction cosine of the normal to the plane containing nodes
+        `n1`, `n2`, and `n3`, reported in the global Cartesian
+        coordinate system.
+
+        Parameters
+        ----------
+        n1 : int
+            Node number
+
+        n2 : int
+            Node number
+
+        n3 : int
+            Node number
+
+        Returns
+        -------
+        float
+            X-direction cosine of the normal
+
+        Examples
+        --------
+        Here we create three nodes in the y-z plane and interrogate the
+        x-component of the normal to that plane, which is trivially 1.0.
+
+        >>> from ansys.mapdl.core import launch_mapdl
+        >>> from ansys.mapdl.core.inline_functions import Query
+        >>> mapdl = launch_mapdl()
+        >>> mapdl.prep7()
+        >>> n1 = mapdl.n(1, 0, 0, 0)
+        >>> n2 = mapdl.n(2, 0, 1, 0)
+        >>> n3 = mapdl.n(3, 0, 1, 1)
+        >>> q = Query(mapdl)
+        >>> q.normnx(n1, n2, n3)
+        1.0
+        """
+        response = self._mapdl.run(f'_=NORMNX({n1}, {n2}, {n3})')
+        float_ = self._parse_parameter_float_response(response)
+        return float_
+
+    def normny(self, n1: int, n2: int, n3: int) -> float:
+        """Y-direction cosine of the normal to the plane containing the given nodes.
+
+        Y-direction cosine of the normal to the plane containing nodes
+        `n1`, `n2`, and `n3`, reported in the global Cartesian
+        coordinate system.
+
+        Parameters
+        ----------
+        n1 : int
+            Node number
+
+        n2 : int
+            Node number
+
+        n3 : int
+            Node number
+
+        Returns
+        -------
+        float
+            Y-direction cosine of the normal
+
+        Examples
+        --------
+        Here we create three nodes in the x-z plane and interrogate the
+        y-component of the normal to that plane, which is trivially 1.0.
+
+        >>> from ansys.mapdl.core import launch_mapdl
+        >>> from ansys.mapdl.core.inline_functions import Query
+        >>> mapdl = launch_mapdl()
+        >>> mapdl.prep7()
+        >>> n1 = mapdl.n(1, 0, 0, 0)
+        >>> n2 = mapdl.n(2, 1, 0, 0)
+        >>> n3 = mapdl.n(3, 1, 0, 1)
+        >>> q = Query(mapdl)
+        >>> q.normny(n1, n2, n3)
+        1.0
+        """
+        response = self._mapdl.run(f'_=NORMNY({n1}, {n2}, {n3})')
+        float_ = self._parse_parameter_float_response(response)
+        return float_
+
+    def normnz(self, n1: int, n2: int, n3: int) -> float:
+        """Z-direction cosine of the normal to the plane containing the given nodes.
+
+        Z-direction cosine of the normal to the plane containing nodes
+        `n1`, `n2`, and `n3`, reported in the global Cartesian
+        coordinate system.
+
+        Parameters
+        ----------
+        n1 : int
+            Node number
+
+        n2 : int
+            Node number
+
+        n3 : int
+            Node number
+
+        Returns
+        -------
+        float
+            Z-direction cosine of the normal
+
+        Examples
+        --------
+        Here we create three nodes in the x-y plane and interrogate the
+        z-component of the normal to that plane, which is trivially 1.0.
+
+        >>> from ansys.mapdl.core import launch_mapdl
+        >>> from ansys.mapdl.core.inline_functions import Query
+        >>> mapdl = launch_mapdl()
+        >>> mapdl.prep7()
+        >>> n1 = mapdl.n(1, 0, 0, 0)
+        >>> n2 = mapdl.n(2, 0, 1, 0)
+        >>> n3 = mapdl.n(3, 1, 1, 0)
+        >>> q = Query(mapdl)
+        >>> q.normnz(n1, n2, n3)
+        1.0
+        """
+        response = self._mapdl.run(f'_=NORMNZ({n1}, {n2}, {n3})')
+        float_ = self._parse_parameter_float_response(response)
+        return float_
+
+
+class _KeypointNormalQueries(_ParameterParsing):
+    _mapdl = None
+
+    def normkx(self, k1: int, k2: int, k3: int) -> float:
+        """X-direction cosine of the normal to the plane containing the given keypoints.
+
+        X-direction cosine of the normal to the plane containing
+        keypoints `k1`, `k2`, and `k3`, reported in the global Cartesian
+        coordinate system.
+
+        Parameters
+        ----------
+        k1 : int
+            Node number
+
+        k2 : int
+            Node number
+
+        k3 : int
+            Node number
+
+        Returns
+        -------
+        float
+            X-direction cosine of the normal
+
+        Examples
+        --------
+        Here we create three keypoints in the y-z plane and interrogate
+        the x-component of the normal to that plane, which is trivially
+        1.0.
+
+        >>> from ansys.mapdl.core import launch_mapdl
+        >>> from ansys.mapdl.core.inline_functions import Query
+        >>> mapdl = launch_mapdl()
+        >>> mapdl.prep7()
+        >>> k1 = mapdl.k(1, 0, 0, 0)
+        >>> k2 = mapdl.k(2, 0, 1, 0)
+        >>> k3 = mapdl.k(3, 0, 1, 1)
+        >>> q = Query(mapdl)
+        >>> q.normnx(k1, k2, k3)
+        1.0
+        """
+        response = self._mapdl.run(f'_=NORMKX({k1}, {k2}, {k3})')
+        float_ = self._parse_parameter_float_response(response)
+        return float_
+
+    def normky(self, k1: int, k2: int, k3: int) -> float:
+        """Y-direction cosine of the normal to the plane containing the given keypoints.
+
+        Y-direction cosine of the normal to the plane containing
+        keypoints `k1`, `k2`, and `k3`, reported in the global Cartesian
+        coordinate system.
+
+        Parameters
+        ----------
+        k1 : int
+            Node number
+
+        k2 : int
+            Node number
+
+        k3 : int
+            Node number
+
+        Returns
+        -------
+        float
+            Y-direction cosine of the normal
+
+        Examples
+        --------
+        Here we create three keypoints in the x-z plane and interrogate
+        the y-component of the normal to that plane, which is trivially
+        1.0.
+
+        >>> from ansys.mapdl.core import launch_mapdl
+        >>> from ansys.mapdl.core.inline_functions import Query
+        >>> mapdl = launch_mapdl()
+        >>> mapdl.prep7()
+        >>> k1 = mapdl.k(1, 0, 0, 0)
+        >>> k2 = mapdl.k(2, 1, 0, 0)
+        >>> k3 = mapdl.k(3, 1, 0, 1)
+        >>> q = Query(mapdl)
+        >>> q.normny(k1, k2, k3)
+        1.0
+        """
+        response = self._mapdl.run(f'_=NORMKY({k1}, {k2}, {k3})')
+        float_ = self._parse_parameter_float_response(response)
+        return float_
+
+    def normkz(self, k1: int, k2: int, k3: int) -> float:
+        """Z-direction cosine of the normal to the plane containing the given keypoints.
+
+        Z-direction cosine of the normal to the plane containing
+        keypoints `k1`, `k2`, and `k3`, reported in the global Cartesian
+        coordinate system.
+
+        Parameters
+        ----------
+        k1 : int
+            Node number
+
+        k2 : int
+            Node number
+
+        k3 : int
+            Node number
+
+        Returns
+        -------
+        float
+            Z-direction cosine of the normal
+
+        Examples
+        --------
+        Here we create three keypoints in the x-y plane and interrogate
+        the z-component of the normal to that plane, which is trivially
+        1.0.
+
+        >>> from ansys.mapdl.core import launch_mapdl
+        >>> from ansys.mapdl.core.inline_functions import Query
+        >>> mapdl = launch_mapdl()
+        >>> mapdl.prep7()
+        >>> k1 = mapdl.k(1, 0, 0, 0)
+        >>> k2 = mapdl.k(2, 0, 1, 0)
+        >>> k3 = mapdl.k(3, 1, 1, 0)
+        >>> q = Query(mapdl)
+        >>> q.normnz(k1, k2, k3)
+        1.0
+        """
+        response = self._mapdl.run(f'_=NORMKZ({k1}, {k2}, {k3})')
+        float_ = self._parse_parameter_float_response(response)
+        return float_
+
diff --git a/ansys/mapdl/core/launcher.py b/ansys/mapdl/core/launcher.py
index 7fe5123f3d..cbaf282b97 100644
--- a/ansys/mapdl/core/launcher.py
+++ b/ansys/mapdl/core/launcher.py
@@ -177,10 +177,10 @@ def launch_grpc(exec_file='', jobname='file', nproc=2, ram=None,
         port available after (or including) this port.
 
     additional_switches : str, optional
-        Additional switches for MAPDL, for example ``"-aa_r"``, and
+        Additional switches for MAPDL, for example ``"-p aa_r"``, the
         academic research license, would be added with:
 
-        - ``additional_switches="-aa_r"``
+        - ``additional_switches="-p aa_r"``
 
         Avoid adding switches like ``"-i"`` ``"-o"`` or ``"-b"`` as
         these are already included to start up the MAPDL server.  See
@@ -208,7 +208,7 @@ def launch_grpc(exec_file='', jobname='file', nproc=2, ram=None,
     --------
     Launch MAPDL using the default configuration.
 
-    >>> from ansys.mapdl import launch_mapdl
+    >>> from ansys.mapdl.core import launch_mapdl
     >>> mapdl = launch_mapdl()
 
     Run MAPDL with shared memory parallel and specify the location of
diff --git a/ansys/mapdl/core/mapdl.py b/ansys/mapdl/core/mapdl.py
index 4b64224719..540e06710b 100644
--- a/ansys/mapdl/core/mapdl.py
+++ b/ansys/mapdl/core/mapdl.py
@@ -19,6 +19,7 @@
 from ansys.mapdl.core.plotting import general_plotter
 from ansys.mapdl.core.post import PostProcessing
 from ansys.mapdl.core.commands import Commands
+from ansys.mapdl.core.inline_functions import Query
 
 
 _PERMITTED_ERRORS = [
@@ -112,6 +113,7 @@ def __init__(self, loglevel='DEBUG', use_vtk=True, log_apdl=False,
                  local=True, **start_parm):
         """Initialize connection with MAPDL. """
         self._show_matplotlib_figures = True  # for testing
+        self._query = None
         self._exited = False
         self._allow_ignore = False
         self._apdl_log = None
@@ -144,6 +146,78 @@ def __init__(self, loglevel='DEBUG', use_vtk=True, log_apdl=False,
 
         self._post = PostProcessing(self)
 
+    def query(self):
+        """Get instance of Query class containing inline functions of APDL.
+
+        Most of the results of these methods are shortcuts for specific
+        combinations of arguments supplied to :func:`ansys.mapdl.core.Mapdl.get`.
+
+        Currently implemented functions:
+
+        - ``centrx(e)`` - get the centroid x-coordinate of element `e`
+        - ``centry(e)`` - get the centroid y-coordinate of element `e`
+        - ``centrz(e)`` - get the centroid z-coordinate of element `e`
+        - ``nx(n)`` - get the x-coordinate of node `n`
+        - ``ny(n)`` - get the y-coordinate of node `n`
+        - ``nz(n)`` - get the z-coordinate of node `n`
+        - ``kx(k)`` - get the x-coordinate of keypoint `k`
+        - ``ky(k)`` - get the y-coordinate of keypoint `k`
+        - ``kz(k)`` - get the z-coordinate of keypoint `k`
+        - ``lx(n, lfrac)`` - X-coordinate of line ``n`` at length fraction ``lfrac``
+        - ``ly(n, lfrac)`` - Y-coordinate of line ``n`` at length fraction ``lfrac``
+        - ``lz(n, lfrac)`` - Z-coordinate of line ``n`` at length fraction ``lfrac``
+        - ``lsx(n, lfrac)`` - X-slope of line ``n`` at length fraction ``lfrac``
+        - ``lsy(n, lfrac)`` - Y-slope of line ``n`` at length fraction ``lfrac``
+        - ``lsz(n, lfrac)`` - Z-slope of line ``n`` at length fraction ``lfrac``
+        - ``ux(n)`` - get the structural displacement at node `n` in x
+        - ``uy(n)`` - get the structural displacement at node `n` in y
+        - ``uz(n)`` - get the structural displacement at node `n` in z
+        - ``rotx(n)`` - get the rotational displacement at node `n` in x
+        - ``roty(n)`` - get the rotational displacement at node `n` in y
+        - ``rotz(n)`` - get the rotational displacement at node `n` in z
+        - ``nsel(n)`` - get the selection status of node `n`
+        - ``ksel(k)`` - get the selection status of keypoint `k`
+        - ``lsel(n)`` - get the selection status of line `n`
+        - ``asel(a)`` - get the selection status of area `a`
+        - ``esel(n)`` - get the selection status of element `e`
+        - ``vsel(v)`` - get the selection status of volume `v`
+        - ``ndnext(n)`` - get the next selected node with a number greater than `n`.
+        - ``kpnext(k)`` - get the next selected keypoint with a number greater than `k`.
+        - ``lsnext(n)`` - get the next selected line with a number greater than `n`.
+        - ``arnext(a)`` - get the next selected area with a number greater than `a`.
+        - ``elnext(e)`` - get the next selected element with a number greater than `e`.
+        - ``vlnext(v)`` - get the next selected volume with a number greater than `v`.
+        - ``node(x, y, z)`` - get the node closest to coordinate (x, y, z)
+        - ``kp(x, y, z)`` - get the keypoint closest to coordinate (x, y, z)
+
+        Returns
+        -------
+        :class:`ansys.mapdl.core.inline_functions.Query`
+            Instance of the Query class
+
+        Examples
+        --------
+        In this example we construct a solid box and mesh it. Then we use
+        the ``Query`` methods ``nx``, ``ny``, and ``nz`` to find the
+        cartesian coordinates of the first node.
+
+        >>> from ansys.mapdl.core import launch_mapdl
+        >>> mapdl = launch_mapdl()
+        >>> mapdl.prep7()
+        >>> mapdl.et(1, 'SOLID5')
+        >>> mapdl.block(0, 10, 0, 20, 0, 30)
+        >>> mapdl.esize(2)
+        >>> mapdl.vmesh('ALL')
+        >>> q = mapdl.query()
+        >>> q.nx(1), q.ny(1), q.nz(1)
+        0.0 20.0 0.0
+
+
+        """
+        if self._query is None:
+            self._query = Query(self)
+        return self._query
+
     @property
     def non_interactive(self):
         """Non-interactive context manager.
@@ -708,7 +782,7 @@ def _close_apdl_log(self):
             self._apdl_log.close()
         self._apdl_log = None
 
-    def nplot(self, knum="", vtk=None, **kwargs):
+    def nplot(self, nnum="", vtk=None, **kwargs):
         """APDL Command: NPLOT
 
         Displays nodes.
@@ -719,12 +793,15 @@ def nplot(self, knum="", vtk=None, **kwargs):
 
         Parameters
         ----------
-        knum : bool, int, optional
+        nnum : bool, int, optional
             Node number key:
 
             - ``False`` : No node numbers on display (default).
             - ``True`` : Include node numbers on display.
 
+            .. note::
+               This parameter is only valid when ``vtk==True``
+
         vtk : bool, optional
             Plot the currently selected nodes using ``pyvista``.
             Defaults to current ``use_vtk`` setting as set on the
@@ -738,7 +815,7 @@ def nplot(self, knum="", vtk=None, **kwargs):
         >>> mapdl.n(1, 0, 0, 0)
         >>> mapdl.n(11, 10, 0, 0)
         >>> mapdl.fill(1, 11, 9)
-        >>> mapdl.nplot(knum=True, vtk=True, background='w', color='k',
+        >>> mapdl.nplot(nnum=True, vtk=True, background='w', color='k',
                         show_bounds=True)
 
         Plot without using VTK
@@ -760,6 +837,9 @@ def nplot(self, knum="", vtk=None, **kwargs):
         if vtk is None:
             vtk = self._use_vtk
 
+        if 'knum' in kwargs:
+            raise ValueError('`knum` keyword deprecated.  Please use `nnum` instead.')
+
         if vtk:
             kwargs.setdefault('title', 'MAPDL Node Plot')
             if not self.mesh.n_node:
@@ -767,7 +847,7 @@ def nplot(self, knum="", vtk=None, **kwargs):
                 return general_plotter([], [], [], **kwargs)
 
             labels = []
-            if knum:
+            if nnum:
                 # must eliminate duplicate points or labeling fails miserably.
                 pcloud = pv.PolyData(self.mesh.nodes)
                 pcloud['labels'] = self.mesh.nnum
@@ -778,11 +858,11 @@ def nplot(self, knum="", vtk=None, **kwargs):
             return general_plotter([], points, labels, **kwargs)
 
         # otherwise, use the built-in nplot
-        if isinstance(knum, bool):
-            knum = int(knum)
+        if isinstance(nnum, bool):
+            nnum = int(nnum)
 
         self._enable_interactive_plotting()
-        return super().nplot(knum, **kwargs)
+        return super().nplot(nnum, **kwargs)
 
     def vplot(self, nv1="", nv2="", ninc="", degen="", scale="",
               vtk=None, quality=4, show_area_numbering=False,
diff --git a/ansys/mapdl/core/math.py b/ansys/mapdl/core/math.py
index 818025cf15..af56d26cdd 100755
--- a/ansys/mapdl/core/math.py
+++ b/ansys/mapdl/core/math.py
@@ -19,7 +19,7 @@
 
 MYCTYPE = {np.int32: 'I',
            np.int64: 'L',
-           np.single: 'S',
+           np.single: 'F',
            np.double: 'D',
            np.complex64: 'C',
            np.complex128: 'Z'}
@@ -59,6 +59,7 @@ def get_nparray_chunks(name, array, chunk_size=DEFAULT_FILE_CHUNK_SIZE):
                                          chunk=chunk)
         i += chunk_size
 
+
 def get_nparray_chunks_mat(name, array, chunk_size=DEFAULT_FILE_CHUNK_SIZE):
     """Serializes a 2D numpy array into chunks
 
@@ -86,21 +87,16 @@ def list_allowed_dtypes():
     return '\n'.join([f'{dtype}' for dtype in dtypes])
 
 
-class MapdlMath():
+class MapdlMath:
     """Abstract mapdl math class.  Created from a ``Mapdl`` instance.
 
     Examples
     --------
-    Create an instance from an existing mapdl instance
-
-    >>> import ansys
-    >>> mapdl = ansys.Mapdl()
-    >>> mm = mapdl.math()
-
-    Alternatively:
+    Create an instance.
 
-    >>> from ansys.mapdl.math import MapdlMath
-    >>> mm = MapdlMath(mapdl)  # from mapdl above
+    >>> from ansys.mapdl.core import launch_mapdl
+    >>> mapdl = launch_mapdl()
+    >>> mm = mapdl.math
 
     Vector addition
 
@@ -108,7 +104,7 @@ class MapdlMath():
     >>> v2 = mm.ones(10)
     >>> v3 = v1 + v2
 
-    Matrix multiplcation
+    Matrix multiplcation (not yet available)
 
     >>> v1 = mm.ones(10)
     >>> m1 = mm.rand(10, 10)
@@ -414,7 +410,7 @@ def matrix(self, matrix, name=None, triu=False):
             return AnsDenseMat(name, self._mapdl)
 
     def load_matrix_from_file(self, dtype=np.double, fname="file.full",
-                              mat_id="STIFF"):
+                              mat_id="STIFF", asarray=False):
         """Import a matrix from an existing FULL file.
 
         Parameters
@@ -440,6 +436,8 @@ def load_matrix_from_file(self, dtype=np.double, fname="file.full",
             * ``"K"``_RE - Real part of the stiffness matrix
             * ``"K"``_IM - Imaginary part of the stiffness matrix
 
+        asarray : bool, optional
+            Return a `scipy` array rather than an APDLMath matrix.
 
         """
         name = id_generator()
@@ -447,9 +445,12 @@ def load_matrix_from_file(self, dtype=np.double, fname="file.full",
                               mat_id, fname)
         self._mapdl.run(f"*SMAT,{name},{MYCTYPE[dtype]},IMPORT,FULL,{fname},{mat_id}",
                         mute=True)
-        return AnsSparseMat(name, self._mapdl)
+        ans_sparse_mat = AnsSparseMat(name, self._mapdl)
+        if asarray:
+            return self._mapdl._mat_data(ans_sparse_mat.id)
+        return ans_sparse_mat
 
-    def stiff(self, dtype=np.double, fname="file.full"):
+    def stiff(self, dtype=np.double, fname="file.full", asarray=False):
         """Load the stiffness matrix from a full file.
 
         Parameters
@@ -460,6 +461,9 @@ def stiff(self, dtype=np.double, fname="file.full"):
         fname : str, optional
             Filename to read the matrix from.
 
+        asarray : bool, optional
+            Return a `scipy` array rather than an APDLMath matrix.
+
         Examples
         --------
         >>> k = mapdl.math.stiff()
@@ -472,9 +476,9 @@ def stiff(self, dtype=np.double, fname="file.full"):
         <60x60 sparse matrix of type '<class 'numpy.float64'>'
             with 1734 stored elements in Compressed Sparse Row format>
         """
-        return self.load_matrix_from_file(dtype, fname, "STIFF")
+        return self.load_matrix_from_file(dtype, fname, "STIFF", asarray)
 
-    def mass(self, dtype=np.double, fname="file.full"):
+    def mass(self, dtype=np.double, fname="file.full", asarray=False):
         """Load the mass matrix from a full file.
 
         Parameters
@@ -485,6 +489,9 @@ def mass(self, dtype=np.double, fname="file.full"):
         fname : str, optional
             Filename to read the matrix from.
 
+        asarray : bool, optional
+            Return a `scipy` array rather than an APDLMath matrix.
+
         Examples
         --------
         >>> m = mapdl.math.mass()
@@ -498,9 +505,9 @@ def mass(self, dtype=np.double, fname="file.full"):
         <60x60 sparse matrix of type '<class 'numpy.float64'>'
             with 1734 stored elements in Compressed Sparse Row format>
         """
-        return self.load_matrix_from_file(dtype, fname, "MASS")
+        return self.load_matrix_from_file(dtype, fname, "MASS", asarray)
 
-    def damp(self, dtype=np.double, fname="file.full"):
+    def damp(self, dtype=np.double, fname="file.full", asarray=False):
         """Load the damping matrix from the full file
 
         Parameters
@@ -511,6 +518,9 @@ def damp(self, dtype=np.double, fname="file.full"):
         fname : str, optional
             Filename to read the matrix from.
 
+        asarray : bool, optional
+            Return a `scipy` array rather than an APDLMath matrix.
+
         Examples
         --------
         >>> m = mapdl.math.damp()
@@ -525,9 +535,9 @@ def damp(self, dtype=np.double, fname="file.full"):
             with 1734 stored elements in Compressed Sparse Row format>
 
         """
-        return self.load_matrix_from_file(dtype, fname, "DAMP")
+        return self.load_matrix_from_file(dtype, fname, "DAMP", asarray)
 
-    def get_vec(self, dtype=np.double, fname="file.full", mat_id="RHS"):
+    def get_vec(self, dtype=np.double, fname="file.full", mat_id="RHS", asarray=False):
         """Load a vector from a file.
 
         Parameters
@@ -547,6 +557,9 @@ def get_vec(self, dtype=np.double, fname="file.full", mat_id="RHS"):
             * ``"BACK"`` - nodal mapping vector (internal to user)
             * ``"FORWARD"`` - nodal mapping vector (user to internal)
 
+        asarray : bool, optional
+            Return a `scipy` array rather than an APDLMath matrix.
+
         Returns
         --------
         ansys.mapdl.math.AnsVec
@@ -564,7 +577,10 @@ def get_vec(self, dtype=np.double, fname="file.full", mat_id="RHS"):
                               mat_id, fname)
         self._mapdl.run(f"*VEC,{name},{MYCTYPE[dtype]},IMPORT,FULL,{fname},{mat_id}",
                         mute=True)
-        return AnsVec(name, self._mapdl)
+        ans_vec = AnsVec(name, self._mapdl)
+        if asarray:
+            return self._mapdl._vec_data(ans_vec.id)
+        return ans_vec
 
     def set_vec(self, data, name=None):
         """Push a numpy array or Python list to the MAPDL Memory
@@ -597,7 +613,7 @@ def set_vec(self, data, name=None):
         self._set_vec(name, data)
         return AnsVec(name, self._mapdl)
 
-    def rhs(self, dtype=np.double, fname="file.full"):
+    def rhs(self, dtype=np.double, fname="file.full", asarray=False):
         """Return the load vector from a full file.
 
         Parameters
@@ -608,6 +624,9 @@ def rhs(self, dtype=np.double, fname="file.full"):
         fname : str, optional
             Filename to read the vector from.  Defaults to ``"file.full"``.
 
+        asarray : bool, optional
+            Return a `scipy` array rather than an APDLMath matrix.
+
         Returns
         -------
         ansys.mapdl.math.AnsVec
@@ -619,7 +638,7 @@ def rhs(self, dtype=np.double, fname="file.full"):
         APDLMath Vector Size 126
 
         """
-        return self.get_vec(dtype, fname, "RHS")
+        return self.get_vec(dtype, fname, "RHS", asarray)
 
     def svd(self, mat, thresh='', sig='', v='', **kwargs):
         """Apply an SVD Algorithm on a matrix.
@@ -703,7 +722,7 @@ def sparse(self, mat, thresh='', **kwargs):
             value is 1E-16.
         """
         kwargs.setdefault('mute', True)
-        self._mapdl.run(f"*COMP,{M.id},SPARSE,{thresh}", **kwargs)
+        self._mapdl.run(f"*COMP,{mat.id},SPARSE,{thresh}", **kwargs)
 
     def eigs(self, nev, k, m=None, c=None, phi=None, algo=None,
              fmin=None, fmax=None):
@@ -997,7 +1016,10 @@ def _send_sparse(self, mname, arr, sym, dtype, chunk_size):
 
         # data vector
         dataname = f'{mname}_DATA'
-        self.set_vec(arr.data, dataname)
+        ans_vec = self.set_vec(arr.data, dataname)
+        if dtype is None:
+            info = self._mapdl._data_info(ans_vec.id)
+            dtype = ANSYS_VALUE_TYPE[info.stype]
 
         # indptr vector
         indptrname = f'{mname}_IND'
@@ -1010,10 +1032,9 @@ def _send_sparse(self, mname, arr, sym, dtype, chunk_size):
         self.set_vec(idx, indxname)
 
         flagsym = 'TRUE' if sym else 'FALSE'
-        self._mapdl.run(f'*SMAT,{mname},D,ALLOC,CSR,{indptrname},{indxname},'
+        self._mapdl.run(f'*SMAT,{mname},{MYCTYPE[dtype]},ALLOC,CSR,{indptrname},{indxname},'
                         f'{dataname},{flagsym}')
 
-
 class ApdlMathObj:
     """Common class for MAPDL Math objects"""
 
@@ -1243,8 +1264,9 @@ def asarray(self):
 
         Examples
         --------
-        >>> from ansys.mapdl import Mapdl
-        >>> mm = mapdl.math()
+        >>> from ansys.mapdl.core import launch_mapdl
+        >>> mapdl = launch_mapdl()
+        >>> mm = mapdl.math
         >>> v = mm.ones(10)
         >>> v.asarray()
         [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
@@ -1290,11 +1312,13 @@ def asarray(self):
 
         Examples
         --------
-        >>> from ansys.mapdl import Mapdl
-        >>> mm = mapdl.math()
+        >>> from ansys.mapdl.core import launch_mapdl
+        >>> mapdl = launch_mapdl()
+        >>> mm = mapdl.math
         >>> v = mm.ones(10)
         >>> v.asarray()
         [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
+        
         """
         return self._mapdl._mat_data(self.id)
 
@@ -1351,8 +1375,9 @@ def T(self):
 
         Examples
         --------
-        >>> from ansys.mapdl import Mapdl
-        >>> mm = mapdl.math()
+        >>> from ansys.mapdl.core import launch_mapdl
+        >>> mapdl = launch_mapdl()
+        >>> mm = mapdl.math
         >>> mat = mm.rand(2, 3)
         >>> mat_t = mat.T
 
diff --git a/ansys/mapdl/core/plotting.py b/ansys/mapdl/core/plotting.py
index cffcff84e9..ec4597baaf 100644
--- a/ansys/mapdl/core/plotting.py
+++ b/ansys/mapdl/core/plotting.py
@@ -15,6 +15,7 @@ def general_plotter(meshes, points, labels,
                     window_size=None,
                     notebook=None,
                     # add_mesh kwargs:
+                    style=None,
                     color='w',
                     show_edges=None, edge_color=None, point_size=5.0,
                     line_width=None, opacity=1.0, flip_scalars=False,
@@ -70,7 +71,13 @@ def general_plotter(meshes, points, labels,
     savefig : str, optional
         Saves screenshot to a file path.
 
-    color : string or 3 item list, optional, defaults to white
+    style : string, optional
+        Visualization style of the mesh.  One of the following:
+        ``style='surface'``, ``style='wireframe'``,
+        ``style='points'``.  Defaults to ``'surface'``. Note that
+        ``'wireframe'`` only shows a wireframe of the outer geometry.
+
+    color : string or 3 item list, optional
         Use to make the entire mesh have a single solid color.  Either
         a string, RGB list, or hex color string.  For example:
         ``color='white'``, ``color='w'``, ``color=[1, 1, 1]``, or
@@ -197,6 +204,7 @@ def general_plotter(meshes, points, labels,
                     scalars=mesh.get('scalars'),
                     scalar_bar_args=scalar_bar_args,
                     color=mesh.get('color', color),
+                    style=mesh.get('style', style),
                     show_edges=show_edges, edge_color=edge_color,
                     smooth_shading=smooth_shading,
                     point_size=point_size, line_width=line_width,
diff --git a/doc/source/user_guide/mesh_geometry.rst b/doc/source/user_guide/mesh_geometry.rst
index 12d555dc22..5f41ebdb70 100644
--- a/doc/source/user_guide/mesh_geometry.rst
+++ b/doc/source/user_guide/mesh_geometry.rst
@@ -69,5 +69,5 @@ commands for creating geometries.
 
 API Reference
 ~~~~~~~~~~~~~
-For a full descrption of the ``Mesh`` and ``Geometry`` classes, please
+For a full description of the ``Mesh`` and ``Geometry`` classes, please
 see :ref:`ref_mesh_api` and :ref:`ref_geometry_api`.
diff --git a/examples/00-mapdl-examples/contact_elements.py b/examples/00-mapdl-examples/contact_elements.py
new file mode 100644
index 0000000000..d7c50faf91
--- /dev/null
+++ b/examples/00-mapdl-examples/contact_elements.py
@@ -0,0 +1,60 @@
+"""
+.. _ref_contact_example:
+
+This example demonstrates how to create contact elements for general
+contact.
+
+Begin by launching MAPDL.
+
+"""
+from ansys.mapdl import core as pymapdl
+
+mapdl = pymapdl.launch_mapdl()
+
+###############################################################################
+# Enter the pre-processor, create a block and mesh it with tetrahedral
+# elements.
+# 
+mapdl.prep7()
+
+vnum0 = mapdl.block(0, 1, 0, 1, 0, 0.5)
+
+mapdl.et(1, 187)
+mapdl.esize(0.1)
+
+mapdl.vmesh(vnum0)
+mapdl.eplot()
+
+###############################################################################
+# Second a volume block above the existing block and mesh it with
+# quadratic hexahedral elements.  Ensure that these blocks do not
+# touch by starting it slightly higher than the existing block.
+#
+# Note how these two blocks do not touch and the mesh is non-conformal.
+
+mapdl.esize(0.09)
+mapdl.et(2, 186)
+mapdl.type(2)
+vnum1 = mapdl.block(0, 1, 0, 1, 0.50001, 1)
+mapdl.vmesh(vnum1)
+mapdl.eplot()
+
+
+###############################################################################
+# Select all the elements at the intersection between the two blocks
+# and generate contact elements.
+
+mapdl.nsel('s', 'loc', 'z', 0.5, 0.50001)
+mapdl.esln('s')
+output = mapdl.gcgen('NEW', splitkey='SPLIT', selopt='SELECT')
+print(output)
+
+###############################################################################
+# Plot the contact element pairs.  Note from the command output above
+# that the section IDs are 5 and 6.
+#
+# Here, we plot the element mesh as a wire-frame to show that the
+# contact pairs overlap.
+
+mapdl.esel('S', 'SEC', vmin=5, vmax=6)
+mapdl.eplot(style='wireframe', line_width=3)
diff --git a/examples/00-mapdl-examples/mapdl_beam.py b/examples/00-mapdl-examples/mapdl_beam.py
index 0514a8fd21..d90284de7a 100644
--- a/examples/00-mapdl-examples/mapdl_beam.py
+++ b/examples/00-mapdl-examples/mapdl_beam.py
@@ -44,7 +44,7 @@
 print(mapdl.mesh.nnum)
 
 # plot the nodes using VTK
-mapdl.nplot(vtk=True, knum=True, cpos='xy', show_bounds=True, point_size=10)
+mapdl.nplot(vtk=True, nnum=True, cpos='xy', show_bounds=True, point_size=10)
 
 ###############################################################################
 # create elements between the nodes
diff --git a/examples/00-mapdl-examples/statically_indeterminate_reaction_force_analysis.py b/examples/00-mapdl-examples/statically_indeterminate_reaction_force_analysis.py
index f4a2139cdc..6ea33121fe 100644
--- a/examples/00-mapdl-examples/statically_indeterminate_reaction_force_analysis.py
+++ b/examples/00-mapdl-examples/statically_indeterminate_reaction_force_analysis.py
@@ -3,7 +3,7 @@
 
 Statically Indeterminate Reaction Force Analysis
 ------------------------------------------------
-Problem Descrption:
+Problem Description:
  - A prismatical bar with built-in ends is loaded axially at two
    intermediate cross sections.  Determine the reactions :math:`R_1`
    and :math:`R_2`.
diff --git a/ignore_words.txt b/ignore_words.txt
index 89d406c8ab..fa444dfad0 100644
--- a/ignore_words.txt
+++ b/ignore_words.txt
@@ -18,3 +18,4 @@ sord
 struc
 emiss
 vise
+sur
diff --git a/requirements_docs.txt b/requirements_docs.txt
index fb46e165ce..ba7436f8c2 100644
--- a/requirements_docs.txt
+++ b/requirements_docs.txt
@@ -12,5 +12,4 @@ sphinx-copybutton
 sphinx-gallery>=0.8.1
 sphinx-autodoc-typehints
 pyansys_sphinx_theme
-pypandoc
-Jinja2==2.11.3  # due to An error happened in rendering the page genindex. with v3
+ansys-mapdl-reader<=0.51.3  # 0.51.4-0.51.5 breaks doc build
diff --git a/requirements_style.txt b/requirements_style.txt
index f249615f64..f851247088 100755
--- a/requirements_style.txt
+++ b/requirements_style.txt
@@ -1,3 +1,3 @@
-codespell==2.0.0
+codespell==2.1.0
 
 
diff --git a/tests/test_inline_functions/test_nearest_queries.py b/tests/test_inline_functions/test_nearest_queries.py
new file mode 100644
index 0000000000..e67d98e7ad
--- /dev/null
+++ b/tests/test_inline_functions/test_nearest_queries.py
@@ -0,0 +1,15 @@
+class TestNearestEntityQueries:
+    def test_nnear(self, box_geometry):
+        q, kps, areas, nodes = box_geometry
+        nearest_node = q.nnear(1)
+        assert nearest_node in nodes
+
+    def test_knear(self, box_geometry):
+        q, kps, areas, nodes = box_geometry
+        nearest_keypoint = q.knear(1)
+        assert nearest_keypoint in kps
+
+    def test_enearn(self, box_geometry):
+        q, kps, areas, nodes = box_geometry
+        nearest_element = q.enearn(1)
+        assert nearest_element > 0
diff --git a/tests/test_inline_functions/test_normals_queries.py b/tests/test_inline_functions/test_normals_queries.py
new file mode 100644
index 0000000000..8c9561383b
--- /dev/null
+++ b/tests/test_inline_functions/test_normals_queries.py
@@ -0,0 +1,64 @@
+
+
+class TestNormalsNodeQueries:
+
+    @staticmethod
+    def build_plane(mapdl, plane: str):
+        n1 = mapdl.n(1, 0, 0, 0)
+        if plane == 'xy':
+            n2 = mapdl.n(2, 0, 1, 0)
+            n3 = mapdl.n(3, 1, 1, 0)
+        elif plane == 'xz':
+            n2 = mapdl.n(2, 1, 0, 0)
+            n3 = mapdl.n(3, 1, 0, 1)
+        elif plane == 'yz':
+            n2 = mapdl.n(2, 0, 1, 0)
+            n3 = mapdl.n(3, 0, 1, 1)
+        return n1, n2, n3
+
+    def test_normnx(self, query):
+        nodes = self.build_plane(query._mapdl, 'yz')
+        cosine = query.normnx(*nodes)
+        assert abs(cosine) == 1.
+
+    def test_normny(self, query):
+        nodes = self.build_plane(query._mapdl, 'xz')
+        cosine = query.normny(*nodes)
+        assert abs(cosine) == 1.
+
+    def test_normnz(self, query):
+        nodes = self.build_plane(query._mapdl, 'xy')
+        cosine = query.normnz(*nodes)
+        assert abs(cosine) == 1.
+
+
+class TestNormalsKeypointsQueries:
+
+    @staticmethod
+    def build_plane(mapdl, plane: str):
+        k1 = mapdl.k(1, 0, 0, 0)
+        if plane == 'xy':
+            k2 = mapdl.k(2, 0, 1, 0)
+            k3 = mapdl.k(3, 1, 1, 0)
+        elif plane == 'xz':
+            k2 = mapdl.k(2, 1, 0, 0)
+            k3 = mapdl.k(3, 1, 0, 1)
+        elif plane == 'yz':
+            k2 = mapdl.k(2, 0, 1, 0)
+            k3 = mapdl.k(3, 0, 1, 1)
+        return k1, k2, k3
+
+    def test_normkx(self, query):
+        keypoints = self.build_plane(query._mapdl, 'yz')
+        cosine = query.normkx(*keypoints)
+        assert abs(cosine) == 1.
+
+    def test_normky(self, query):
+        keypoints = self.build_plane(query._mapdl, 'xz')
+        cosine = query.normky(*keypoints)
+        assert abs(cosine) == 1.
+
+    def test_normkz(self, query):
+        keypoints = self.build_plane(query._mapdl, 'xy')
+        cosine = query.normkz(*keypoints)
+        assert abs(cosine) == 1.
diff --git a/tests/test_mapdl.py b/tests/test_mapdl.py
index d32c3b5d51..46005fb954 100644
--- a/tests/test_mapdl.py
+++ b/tests/test_mapdl.py
@@ -360,13 +360,13 @@ def test_enum(mapdl, make_block):
     assert np.allclose(mapdl.mesh.enum, range(1, mapdl.mesh.n_elem + 1))
 
 
-@pytest.mark.parametrize('knum', [True, False])
+@pytest.mark.parametrize('nnum', [True, False])
 @skip_no_xserver
-def test_nplot_vtk(cleared, mapdl, knum):
+def test_nplot_vtk(cleared, mapdl, nnum):
     mapdl.n(1, 0, 0, 0)
     mapdl.n(11, 10, 0, 0)
     mapdl.fill(1, 11, 9)
-    mapdl.nplot(vtk=True, knum=knum, background='w', color='k')
+    mapdl.nplot(vtk=True, nnum=nnum, background='w', color='k')
 
 
 @skip_no_xserver