diff --git a/src/ansys/mapdl/core/mesh/mesh.py b/src/ansys/mapdl/core/mesh/mesh.py
index a4efe3f6cf..03ccef42bb 100644
--- a/src/ansys/mapdl/core/mesh/mesh.py
+++ b/src/ansys/mapdl/core/mesh/mesh.py
@@ -133,12 +133,12 @@ def _parse_vtk(
                 tshape_label = SHAPE_MAP[tshape_num]
                 type_ref[etype_ind] = TARGE170_MAP.get(tshape_label, 0)
 
+    nodes, angles, nnum = mesh.nodes, mesh.node_angles, mesh.nnum
+
     offset, celltypes, cells = _reader.ans_vtk_convert(
-        mesh._elem, mesh._elem_off, type_ref, mesh.nnum, True
+        mesh._elem, mesh._elem_off, type_ref, nnum, True
     )  # for reset_midside
 
-    nodes, angles, nnum = mesh.nodes, mesh.node_angles, mesh.nnum
-
     # fix missing midside
     if np.any(cells == -1):
         if fix_midside:
@@ -158,7 +158,13 @@ def _parse_vtk(
         grid = pv.UnstructuredGrid(offset, cells, celltypes, nodes, deep=True)
 
     # Store original ANSYS element and node information
-    grid.point_data["ansys_node_num"] = nnum
+    try:
+        grid.point_data["ansys_node_num"] = nnum
+    except ValueError:
+        grid.point_data["ansys_node_num"] = (
+            mesh._mapdl.nlist(kinternal="internal").to_array()[:, 0].astype(np.int32)
+        )
+
     grid.cell_data["ansys_elem_num"] = mesh.enum
     grid.cell_data["ansys_real_constant"] = mesh.elem_real_constant
     grid.cell_data["ansys_material_type"] = mesh.material_type
diff --git a/tests/conftest.py b/tests/conftest.py
index 1050bfb0c5..29ec833b46 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -869,3 +869,82 @@ def contact_solve(mapdl):
     mapdl.allsel()
     mapdl.set("last")
     mapdl.mute = False
+
+
+@pytest.fixture(scope="function")
+def cuadratic_beam_problem(mapdl):
+    mapdl.clear()
+
+    # Enter verification example mode and the pre-processing routine.
+    mapdl.prep7()
+
+    # Type of analysis: static.
+    mapdl.antype("STATIC")
+
+    # Element type: BEAM188.
+    mapdl.et(1, "BEAM188")
+
+    # Special Features are defined by keyoptions of beam element:
+
+    # KEYOPT(3)
+    # Shape functions along the length:
+    # Cubic
+    mapdl.keyopt(1, 3, 3)  # Cubic shape function
+
+    # KEYOPT(9)
+    # Output control for values extrapolated to the element
+    # and section nodes:
+    # Same as KEYOPT(9) = 1 plus stresses and strains at all section nodes
+    mapdl.keyopt(1, 9, 3, mute=True)
+
+    mapdl.mp("EX", 1, 30e6)
+    mapdl.mp("PRXY", 1, 0.3)
+    print(mapdl.mplist())
+
+    w_f = 1.048394965
+    w_w = 0.6856481
+    sec_num = 1
+    mapdl.sectype(sec_num, "BEAM", "I", "ISection")
+    mapdl.secdata(15, 15, 28 + (2 * w_f), w_f, w_f, w_w)
+
+    # Define nodes
+    for node_num in range(1, 6):
+        mapdl.n(node_num, (node_num - 1) * 120, 0, 0)
+
+    # Define one node for the orientation of the beam cross-section.
+    orient_node = mapdl.n(6, 60, 1)
+
+    # Print the list of the created nodes.
+    print(mapdl.nlist())
+
+    for elem_num in range(1, 5):
+        mapdl.e(elem_num, elem_num + 1, orient_node)
+
+    # Print the list of the created elements.
+    print(mapdl.elist())
+
+    # Display elements with their nodes numbers.
+    mapdl.eplot(show_node_numbering=True, line_width=5, cpos="xy", font_size=40)
+
+    # BC for the beams seats
+    mapdl.d(2, "UX", lab2="UY")
+    mapdl.d(4, "UY")
+
+    # BC for all nodes of the beam
+    mapdl.nsel("S", "LOC", "Y", 0)
+    mapdl.d("ALL", "UZ")
+    mapdl.d("ALL", "ROTX")
+    mapdl.d("ALL", "ROTY")
+    mapdl.nsel("ALL")
+
+    # Parametrization of the distributed load.
+    w = 10000 / 12
+
+    # Application of the surface load to the beam element.
+    mapdl.sfbeam(1, 1, "PRES", w)
+    mapdl.sfbeam(4, 1, "PRES", w)
+    mapdl.finish()
+
+    mapdl.run("/SOLU")
+    mapdl.solve()
+    mapdl.finish()
diff --git a/tests/test_mapdl.py b/tests/test_mapdl.py
index b0b24255f9..817af3934d 100644
--- a/tests/test_mapdl.py
+++ b/tests/test_mapdl.py
@@ -1890,3 +1890,14 @@ def test_igesin_whitespace(mapdl, cleared, tmpdir):
     out = mapdl.igesin(dest)
     n_ent = re.findall(r"TOTAL NUMBER OF ENTITIES \s*=\s*(\d*)", out)
     assert int(n_ent[0]) > 0
+
+
+def test_cuadratic_beam(mapdl, cuadratic_beam_problem):
+    mapdl.post1()
+    mapdl.set(1)
+    assert (
+        mapdl.post_processing.plot_nodal_displacement(
+            "NORM", line_width=10, render_lines_as_tubes=True, smooth_shading=True
+        )
+        is None
+    )