diff --git a/CHANGELOG b/CHANGELOG
index a328b3f8..6bdf0c0f 100644
--- a/CHANGELOG
+++ b/CHANGELOG
@@ -12,6 +12,13 @@
 Version 4.2.1 (development)
 ===========================
 
+- The GLVis auto refinement algorithm now takes into account the order of the
+  data (mesh and grid function). The new refinement is chosen to be sufficient
+  for resolving the curvature of the data, but only if we can do that with less
+  than 2M vertices and 16 refinements. Otherwise, we print a warning and the
+  user may still need to press 'o' to fully resolve the data. For more details,
+  see the section Auto-refinement in README.md.
+
 - Significantly improved memory usage.
 
 - Add support to visualize solutions on 1D elements embedded in 2D and 3D.
diff --git a/README.md b/README.md
index 530fbe73..ad287e42 100644
--- a/README.md
+++ b/README.md
@@ -68,6 +68,12 @@ is a partial list of the available functionality. Some of these keys can also be
 provided as input, using the `-k` command-line option and the `keys` script
 command.
 
+For high-order meshes and/or solution data, GLVis performs element subdivision
+to try to represent the data more accurately. However, for highly-varying data
+or large meshes, the auto-selected subdivision factor (see the
+[Auto-refinement](#auto-refinement) section below) may not be sufficient -- use
+the keys <kbd>o</kbd> / <kbd>O</kbd>, described below, to manually adjust the
+subdivision factor.
 
 SPDX-License-Identifier: BSD-3-Clause <br>
 LLNL Release Number: LLNL-CODE-443271 <br>
@@ -298,3 +304,46 @@ Key commands
   - `x`-component: `v_x`
   - `y`-component: `v_y`
   - `z`-component: `v_z`
+
+## Auto-refinement
+
+The GLVis auto-refinement algorithm selects a subdivision factor trying to
+achieve an accurate representation of high-order meshes and solution data while
+keeping the initial time to visualize the data reasonable. The algorithm can be
+summarized as follows:
+- GLVis draws surface elements; the number of drawn elements, `ne`, is either:
+  - the number of elements in the mesh for 2D meshes (including surface meshes,
+    i.e. 2D meshes embedded in 3D space), or
+  - the number of boundary mesh elements described by the mesh in 3D.
+- A tentative upper limit on the number of vertices to be drawn is defined based
+  on the maximum order of the mesh and the solution, `max_order`:
+  ```
+  max_vert = ne * (max_order + 1) * (max_order + 1)
+  ```
+- To allow more accurate representation for small meshes, this number is
+  potentially increased:
+  ```
+  max_vert = max(max_vert, 100 000)
+  ```
+- To keep the time to initially visualize the data reasonable, this number is
+  potentially reduced:
+  ```
+  max_vert = min(max_vert, 2 000 000)
+  ```
+- Finally, the subdivision factor `ref` is chosen to be the largest number such
+  that:
+  - the number of vertices needed to draw the `ne` surface elements with `ref`
+    subdivisions does not exceed `max_vert`:
+    ```
+    ne * (ref + 1) * (ref + 1) <= max_vert
+    ```
+  - for large meshes where the above limit cannot be satisfied, set `ref = 1`
+  - for small meshes, avoid excessive refinements:
+    ```
+    ref <= 16
+    ```
+
+Note that, for highly-varying data or large meshes, this auto-selected
+subdivision factor may not be sufficient for accurate representation. In such
+cases the subdivision can be manually adjusted using the keys <kbd>o</kbd> /
+<kbd>O</kbd>, described above.
diff --git a/lib/vsdata.cpp b/lib/vsdata.cpp
index ff3cb835..37338acd 100644
--- a/lib/vsdata.cpp
+++ b/lib/vsdata.cpp
@@ -47,6 +47,60 @@ void VisualizationSceneScalarData::FixValueRange()
    }
 }
 
+int VisualizationSceneScalarData::GetFunctionAutoRefineFactor(GridFunction &gf)
+{
+   Mesh *mesh = gf.FESpace()->GetMesh();
+   const int order = gf.FESpace()->GetMaxElementOrder();
+
+   // check for integral elements
+   const int dim = mesh->Dimension();
+   const FiniteElementCollection *fec = gf.FESpace()->FEColl();
+   if (fec && fec->GetMapType(dim) == FiniteElement::INTEGRAL)
+   {
+      cout << "Warning: integral elements are non-polynomial in the physical space,\n"
+           << "         consider increasing the refinement by the key 'o'."
+           << endl;
+   }
+
+   return std::max(order, 1);
+}
+
+int VisualizationSceneScalarData::GetAutoRefineFactor()
+{
+   const int dim = mesh->Dimension();
+   const int ne = (dim == 3)?(mesh->GetNBE()):(mesh->GetNE());
+
+   // determine the refinement based on the order of the mesh and grid function
+   int order_ref = GetFunctionAutoRefineFactor();
+
+   // mesh
+   const FiniteElementSpace *nfes = mesh->GetNodalFESpace();
+   if (nfes)
+   {
+      const int order = nfes->GetMaxElementOrder();
+      order_ref = std::max(order_ref, order);
+   }
+
+   // limit the total number of vertices
+   int auto_ref_surf_vert = ne * (order_ref+1) * (order_ref+1);
+   auto_ref_surf_vert = std::min(std::max(auto_ref_surf_vert,
+                                          auto_ref_min_surf_vert), auto_ref_max_surf_vert);
+
+   // approach the given number of vertices
+   int ref = 1;
+   while (ref < auto_ref_max && ne*(ref+2)*(ref+2) <= auto_ref_surf_vert)
+   { ref++; }
+
+   if (ref < order_ref)
+   {
+      cout << "Warning: the automatic refinement does not resolve the data fully,\n"
+           << "         consider increasing the refinement by the key 'o'."
+           << endl;
+   }
+
+   return ref;
+}
+
 void VisualizationSceneScalarData::DoAutoscale(bool prepare)
 {
    if (autoscale == 1)
@@ -1321,7 +1375,8 @@ void VisualizationSceneScalarData::Init()
    scaling = 0;
    drawaxes = colorbar = 0;
    auto_ref_max = 16;
-   auto_ref_max_surf_elem = 20000;
+   auto_ref_min_surf_vert = 100000;
+   auto_ref_max_surf_vert = 2000000;
    minv = 0.0;
    maxv = 1.0;
    logscale = false;
diff --git a/lib/vsdata.hpp b/lib/vsdata.hpp
index a9605d43..9d444a30 100644
--- a/lib/vsdata.hpp
+++ b/lib/vsdata.hpp
@@ -78,7 +78,7 @@ class VisualizationSceneScalarData : public VisualizationScene
 
    int scaling, colorbar, drawaxes;
    Shading shading;
-   int auto_ref_max, auto_ref_max_surf_elem;
+   int auto_ref_max, auto_ref_min_surf_vert, auto_ref_max_surf_vert;
 
    // Formatter for axes & colorbar numbers. Set defaults.
    function<string(double)> axis_formatter = NumberFormatter(4, 'd', false);
@@ -143,6 +143,10 @@ class VisualizationSceneScalarData : public VisualizationScene
 
    void FixValueRange();
 
+   static int GetFunctionAutoRefineFactor(GridFunction &gf);
+   virtual int GetFunctionAutoRefineFactor() = 0;
+   virtual int GetAutoRefineFactor();
+
    void Cone(gl3::GlDrawable& buf, glm::mat4 transform, double cval);
 
 public:
@@ -202,13 +206,13 @@ class VisualizationSceneScalarData : public VisualizationScene
    void SetValueRange(double, double);
 
    virtual void SetShading(Shading, bool) = 0;
-   virtual void ToogleShading() { SetShading((Shading)(((int)shading + 1) % (int)Shading::Max), true); }
+   virtual void ToggleShading() { SetShading((Shading)(((int)shading + 1) % (int)Shading::Max), true); }
    virtual Shading GetShading() { return shading; }
    virtual void SetRefineFactors(int, int) = 0;
-   void SetAutoRefineLimits(int max_ref, int max_surf_elem)
+   void SetAutoRefineLimits(int max_ref, int max_surf_vert)
    {
       auto_ref_max = max_ref;
-      auto_ref_max_surf_elem = max_surf_elem;
+      auto_ref_max_surf_vert = max_surf_vert;
    }
    virtual void AutoRefine() = 0;
    virtual void ToggleAttributes(Array<int> &attr_list) = 0;
diff --git a/lib/vssolution.cpp b/lib/vssolution.cpp
index 0e4473da..a529d47f 100644
--- a/lib/vssolution.cpp
+++ b/lib/vssolution.cpp
@@ -589,10 +589,15 @@ void VisualizationSceneSolution::ToggleDrawElems()
 void VisualizationSceneSolution::NewMeshAndSolution(
    Mesh *new_m, Mesh *new_mc, Vector *new_sol, GridFunction *new_u)
 {
+   Mesh *old_m = mesh;
+   mesh = new_m;
+   mesh_coarse = new_mc;
+   sol = new_sol;
+   rsol = new_u;
+
    // If the number of elements changes, recompute the refinement factor
-   if (mesh->GetNE() != new_m->GetNE())
+   if (mesh->GetNE() != old_m->GetNE())
    {
-      mesh = new_m;
       int ref = GetAutoRefineFactor();
       if (TimesToRefine != ref || EdgeRefineFactor != 1)
       {
@@ -601,10 +606,6 @@ void VisualizationSceneSolution::NewMeshAndSolution(
          cout << "Subdivision factors = " << TimesToRefine << ", 1" << endl;
       }
    }
-   mesh = new_m;
-   mesh_coarse = new_mc;
-   sol = new_sol;
-   rsol = new_u;
 
    have_sol_range = false;
    DoAutoscale(false);
@@ -881,7 +882,7 @@ void VisualizationSceneSolution::ToggleShading()
 {
    if (rsol)
    {
-      VisualizationSceneScalarData::ToogleShading();
+      VisualizationSceneScalarData::ToggleShading();
    }
    else
    {
@@ -988,16 +989,11 @@ void VisualizationSceneSolution::SetRefineFactors(int tot, int bdr)
    }
 }
 
-int VisualizationSceneSolution::GetAutoRefineFactor()
+int VisualizationSceneSolution::GetFunctionAutoRefineFactor()
 {
-   int ne = mesh->GetNE(), ref = 1;
-
-   while (ref < auto_ref_max && ne*(ref+1)*(ref+1) <= auto_ref_max_surf_elem)
-   {
-      ref++;
-   }
+   if (!rsol) { return 1; }
 
-   return ref;
+   return VisualizationSceneScalarData::GetFunctionAutoRefineFactor(*rsol);
 }
 
 void VisualizationSceneSolution::AutoRefine()
diff --git a/lib/vssolution.hpp b/lib/vssolution.hpp
index 2d70517a..1e12d136 100644
--- a/lib/vssolution.hpp
+++ b/lib/vssolution.hpp
@@ -77,7 +77,7 @@ class VisualizationSceneSolution : public VisualizationSceneScalarData
                         Vector &values, int sides, Array<double> &lvl,
                         int flat = 0);
 
-   int GetAutoRefineFactor();
+   int GetFunctionAutoRefineFactor() override;
 
    // Used for drawing markers for element and vertex numbering
    double GetElementLengthScale(int k);
@@ -92,35 +92,35 @@ class VisualizationSceneSolution : public VisualizationSceneScalarData
 
    virtual ~VisualizationSceneSolution();
 
-   virtual std::string GetHelpString() const;
+   std::string GetHelpString() const override;
 
    void SetGridFunction(GridFunction & u) { rsol = &u; }
 
    void NewMeshAndSolution(Mesh *new_m, Mesh *new_mc, Vector *new_sol,
                            GridFunction *new_u = NULL);
 
-   virtual void SetNewScalingFromBox();
-   virtual void FindNewBox(bool prepare);
-   virtual void FindNewValueRange(bool prepare);
-   virtual void FindNewBoxAndValueRange(bool prepare)
+   void SetNewScalingFromBox() override;
+   void FindNewBox(bool prepare) override;
+   void FindNewValueRange(bool prepare) override;
+   void FindNewBoxAndValueRange(bool prepare) override
    { FindNewBox(prepare); }
-   virtual void FindMeshBox(bool prepare);
+   void FindMeshBox(bool prepare) override;
 
-   virtual void ToggleLogscale(bool print);
-   virtual void EventUpdateBackground();
-   virtual void EventUpdateColors();
-   virtual void UpdateLevelLines() { PrepareLevelCurves(); }
-   virtual void UpdateValueRange(bool prepare);
+   void ToggleLogscale(bool print) override;
+   void EventUpdateBackground() override;
+   void EventUpdateColors() override;
+   void UpdateLevelLines()  override { PrepareLevelCurves(); }
+   void UpdateValueRange(bool prepare) override;
 
    void PrepareWithNormals();
    void PrepareFlat();
    void PrepareFlat2();
 
-   virtual void PrepareLines();
+   void PrepareLines() override;
    void PrepareLines2();
    void PrepareLines3();
 
-   virtual void Prepare();
+   void Prepare() override;
    void PrepareLevelCurves();
    void PrepareLevelCurves2();
 
@@ -140,12 +140,12 @@ class VisualizationSceneSolution : public VisualizationSceneScalarData
 
    void PrepareCP();
 
-   virtual gl3::SceneInfo GetSceneObjs();
+   gl3::SceneInfo GetSceneObjs() override;
 
    void glTF_ExportBoundary(glTF_Builder &bld,
                             glTF_Builder::buffer_id buffer,
                             glTF_Builder::material_id black_mat);
-   virtual void glTF_Export();
+   void glTF_Export() override;
 
    void ToggleDrawBdr();
 
@@ -164,17 +164,17 @@ class VisualizationSceneSolution : public VisualizationSceneScalarData
       PrepareNumbering(false);
    }
 
-   virtual void SetShading(Shading, bool);
-   virtual void ToggleShading();
+   void SetShading(Shading, bool) override;
+   void ToggleShading() override;
 
    void ToggleDrawCP() { draw_cp = !draw_cp; PrepareCP(); }
 
    void ToggleRefinements();
    void ToggleRefinementFunction();
 
-   virtual void SetRefineFactors(int, int);
-   virtual void AutoRefine();
-   virtual void ToggleAttributes(Array<int> &attr_list);
+   void SetRefineFactors(int, int) override;
+   void AutoRefine() override;
+   void ToggleAttributes(Array<int> &attr_list) override;
 
    virtual void SetDrawMesh(int i) { drawmesh = i % 3; }
    virtual int GetDrawMesh() { return drawmesh; }
diff --git a/lib/vssolution3d.cpp b/lib/vssolution3d.cpp
index 49b74ee0..20452f3a 100644
--- a/lib/vssolution3d.cpp
+++ b/lib/vssolution3d.cpp
@@ -829,12 +829,18 @@ void VisualizationSceneSolution3d::NewMeshAndSolution(
       delete [] node_pos;
       node_pos = new double[new_m->GetNV()];
    }
+
+   Mesh *old_m = mesh;
+   mesh = new_m;
+   mesh_coarse = new_mc;
+   sol = new_sol;
+   GridF = new_u;
+
    // If the number of surface elements changes, recompute the refinement factor
-   if (mesh->Dimension() != new_m->Dimension() ||
-       (mesh->Dimension() == 2 && mesh->GetNE() != new_m->GetNE()) ||
-       (mesh->Dimension() == 3 && mesh->GetNBE() != new_m->GetNBE()))
+   if (mesh->Dimension() != old_m->Dimension() ||
+       (mesh->Dimension() == 2 && mesh->GetNE() != old_m->GetNE()) ||
+       (mesh->Dimension() == 3 && mesh->GetNBE() != old_m->GetNBE()))
    {
-      mesh = new_m;
       int ref = GetAutoRefineFactor();
       if (TimesToRefine != ref)
       {
@@ -842,10 +848,6 @@ void VisualizationSceneSolution3d::NewMeshAndSolution(
          cout << "Subdivision factor = " << TimesToRefine << endl;
       }
    }
-   mesh = new_m;
-   mesh_coarse = new_mc;
-   sol = new_sol;
-   GridF = new_u;
    FindNodePos();
 
    DoAutoscale(false);
@@ -892,7 +894,7 @@ void VisualizationSceneSolution3d::ToggleShading()
 {
    if (GridF)
    {
-      VisualizationSceneScalarData::ToogleShading();
+      VisualizationSceneScalarData::ToggleShading();
    }
    else
    {
@@ -919,20 +921,11 @@ void VisualizationSceneSolution3d::SetRefineFactors(int f, int ignored)
    }
 }
 
-int VisualizationSceneSolution3d::GetAutoRefineFactor()
+int VisualizationSceneSolution3d::GetFunctionAutoRefineFactor()
 {
-   int ne = mesh->GetNBE(), ref = 1;
-   if (mesh->Dimension() == 2)
-   {
-      ne = mesh->GetNE();
-   }
-
-   while (ref < auto_ref_max && ne*(ref+1)*(ref+1) <= auto_ref_max_surf_elem)
-   {
-      ref++;
-   }
+   if (!GridF) { return 1; }
 
-   return ref;
+   return VisualizationSceneScalarData::GetFunctionAutoRefineFactor(*GridF);
 }
 
 void VisualizationSceneSolution3d::AutoRefine()
diff --git a/lib/vssolution3d.hpp b/lib/vssolution3d.hpp
index 53c91a24..ce12bdf1 100644
--- a/lib/vssolution3d.hpp
+++ b/lib/vssolution3d.hpp
@@ -101,7 +101,7 @@ class VisualizationSceneSolution3d : public VisualizationSceneScalarData
                                 const int nh, const int face_splits,
                                 const DenseMatrix *grad = NULL);
 
-   int GetAutoRefineFactor();
+   int GetFunctionAutoRefineFactor() override;
 
    bool CheckPositions(Array<int> &vertices) const
    {
@@ -129,22 +129,22 @@ class VisualizationSceneSolution3d : public VisualizationSceneScalarData
 
    virtual ~VisualizationSceneSolution3d();
 
-   virtual std::string GetHelpString() const;
+   std::string GetHelpString() const override;
 
-   virtual void FindNewBox(bool prepare);
-   virtual void FindNewValueRange(bool prepare);
+   void FindNewBox(bool prepare) override;
+   void FindNewValueRange(bool prepare) override;
 
-   virtual void PrepareRuler()
+   void PrepareRuler() override
    { VisualizationSceneScalarData::PrepareRuler(false); }
    virtual void PrepareFlat();
-   virtual void PrepareLines();
-   virtual void Prepare();
+   void PrepareLines() override;
+   void Prepare() override;
    virtual void PrepareOrderingCurve();
    virtual void PrepareOrderingCurve1(gl3::GlDrawable& buf, bool arrows,
                                       bool color);
-   virtual gl3::SceneInfo GetSceneObjs();
+   gl3::SceneInfo GetSceneObjs() override;
 
-   virtual void glTF_Export();
+   void glTF_Export() override;
 
    void ToggleDrawElems()
    { drawelems = !drawelems; Prepare(); }
@@ -155,11 +155,11 @@ class VisualizationSceneSolution3d : public VisualizationSceneScalarData
    //           3 - no arrows (black), 4 - with arrows (black)
    void ToggleDrawOrdering() { draworder = (draworder+1)%5; }
 
-   virtual void SetShading(Shading, bool);
-   virtual void ToggleShading();
-   virtual void SetRefineFactors(int, int);
-   virtual void AutoRefine();
-   virtual void ToggleAttributes(Array<int> &attr_list);
+   void SetShading(Shading, bool) override;
+   void ToggleShading() override;
+   void SetRefineFactors(int, int) override;
+   void AutoRefine() override;
+   void ToggleAttributes(Array<int> &attr_list) override;
 
    void FindNodePos();
 
@@ -188,10 +188,10 @@ class VisualizationSceneSolution3d : public VisualizationSceneScalarData
    void ToggleCPAlgorithm();
    void MoveLevelSurf(int);
    void NumberOfLevelSurf(int);
-   virtual void EventUpdateColors();
-   virtual void UpdateLevelLines()
+   void EventUpdateColors() override;
+   void UpdateLevelLines() override
    { PrepareLines(); PrepareCuttingPlaneLines(); }
-   virtual void UpdateValueRange(bool prepare);
+   void UpdateValueRange(bool prepare) override;
 
    virtual void SetDrawMesh(int i)
    {
diff --git a/lib/vsvector.cpp b/lib/vsvector.cpp
index 5a3a997e..2e5a3fe8 100644
--- a/lib/vsvector.cpp
+++ b/lib/vsvector.cpp
@@ -418,13 +418,15 @@ void VisualizationSceneVector::NewMeshAndSolution(GridFunction &vgf, Mesh *mc)
       delete solx;
    }
 
+   Mesh *old_m = mesh;
+   Mesh *new_mesh = vgf.FESpace()->GetMesh();
+   mesh = new_mesh;
+   mesh_coarse = mc;
    VecGridF = &vgf;
 
    // If the number of elements changes, recompute the refinement factor
-   Mesh *new_mesh = vgf.FESpace()->GetMesh();
-   if (mesh->GetNE() != new_mesh->GetNE())
+   if (mesh->GetNE() != old_m->GetNE())
    {
-      mesh = new_mesh;
       int ref = GetAutoRefineFactor();
       if (TimesToRefine != ref || EdgeRefineFactor != 1)
       {
@@ -438,8 +440,6 @@ void VisualizationSceneVector::NewMeshAndSolution(GridFunction &vgf, Mesh *mc)
          cout << "Vector subdivision factor = 1" << endl;
       }
    }
-   mesh = new_mesh;
-   mesh_coarse = mc;
 
    solx = new Vector(mesh->GetNV());
    soly = new Vector(mesh->GetNV());
@@ -885,6 +885,13 @@ void VisualizationSceneVector::DrawVector(double px, double py, double vx,
    }
 }
 
+int VisualizationSceneVector::GetFunctionAutoRefineFactor()
+{
+   if (!VecGridF) { return 1;}
+
+   return VisualizationSceneScalarData::GetFunctionAutoRefineFactor(*VecGridF);
+}
+
 void VisualizationSceneVector::PrepareVectorField()
 {
    int rerun;
diff --git a/lib/vsvector.hpp b/lib/vsvector.hpp
index 9ce9707c..46dd8209 100644
--- a/lib/vsvector.hpp
+++ b/lib/vsvector.hpp
@@ -30,11 +30,11 @@ class VisualizationSceneVector : public VisualizationSceneSolution
 
    void Init();
 
-   virtual void GetRefinedValues(int i, const IntegrationRule &ir,
-                                 Vector &vals, DenseMatrix &tr);
-   virtual int GetRefinedValuesAndNormals(int i, const IntegrationRule &ir,
-                                          Vector &vals, DenseMatrix &tr,
-                                          DenseMatrix &normals);
+   void GetRefinedValues(int i, const IntegrationRule &ir,
+                         Vector &vals, DenseMatrix &tr) override;
+   int GetRefinedValuesAndNormals(int i, const IntegrationRule &ir,
+                                  Vector &vals, DenseMatrix &tr,
+                                  DenseMatrix &normals) override;
 
    double (*Vec2Scalar)(double, double);
 
@@ -45,6 +45,8 @@ class VisualizationSceneVector : public VisualizationSceneSolution
    Vector vc0;
    IsoparametricTransformation T0;
 
+   int GetFunctionAutoRefineFactor() override;
+
 public:
    VisualizationSceneVector(Mesh &m, Vector &sx, Vector &sy, Mesh *mc = NULL);
    VisualizationSceneVector(GridFunction &vgf);
@@ -53,14 +55,14 @@ class VisualizationSceneVector : public VisualizationSceneSolution
 
    virtual ~VisualizationSceneVector();
 
-   virtual std::string GetHelpString() const;
+   std::string GetHelpString() const override;
 
    void NPressed();
    void PrepareDisplacedMesh();
-   virtual void PrepareLines()
+   void PrepareLines() override
    { VisualizationSceneSolution::PrepareLines(); PrepareDisplacedMesh(); }
 
-   virtual void ToggleDrawElems();
+   void ToggleDrawElems() override;
 
    virtual void PrepareVectorField();
    void ToggleVectorField();
@@ -74,11 +76,11 @@ class VisualizationSceneVector : public VisualizationSceneSolution
       }
    }
 
-   virtual gl3::SceneInfo GetSceneObjs();
+   gl3::SceneInfo GetSceneObjs() override;
 
-   virtual void glTF_Export();
+   void glTF_Export() override;
 
-   virtual void EventUpdateColors() { Prepare(); PrepareVectorField(); }
+   void EventUpdateColors() override { Prepare(); PrepareVectorField(); }
 
    // refinement factor for the vectors
    int RefineFactor;
diff --git a/lib/vsvector3d.cpp b/lib/vsvector3d.cpp
index abb49bea..15bf0119 100644
--- a/lib/vsvector3d.cpp
+++ b/lib/vsvector3d.cpp
@@ -442,6 +442,13 @@ void VisualizationSceneVector3d::Init()
    }
 }
 
+int VisualizationSceneVector3d::GetFunctionAutoRefineFactor()
+{
+   if (!VecGridF) { return 1; }
+
+   return VisualizationSceneScalarData::GetFunctionAutoRefineFactor(*VecGridF);
+}
+
 VisualizationSceneVector3d::~VisualizationSceneVector3d()
 {
    delete sol;
@@ -474,12 +481,16 @@ void VisualizationSceneVector3d::NewMeshAndSolution(
       node_pos = new double[new_m->GetNV()];
    }
 
+   Mesh *old_m = mesh;
+   VecGridF = new_v;
+   mesh = new_m;
+   mesh_coarse = new_mc;
+
    // If the number of surface elements changes, recompute the refinement factor
-   if (mesh->Dimension() != new_m->Dimension() ||
-       (mesh->Dimension() == 2 && mesh->GetNE() != new_m->GetNE()) ||
-       (mesh->Dimension() == 3 && mesh->GetNBE() != new_m->GetNBE()))
+   if (mesh->Dimension() != old_m->Dimension() ||
+       (mesh->Dimension() == 2 && mesh->GetNE() != old_m->GetNE()) ||
+       (mesh->Dimension() == 3 && mesh->GetNBE() != old_m->GetNBE()))
    {
-      mesh = new_m;
       int ref = GetAutoRefineFactor();
       if (TimesToRefine != ref)
       {
@@ -490,9 +501,6 @@ void VisualizationSceneVector3d::NewMeshAndSolution(
 
    FiniteElementSpace *new_fes = new_v->FESpace();
 
-   VecGridF = new_v;
-   mesh = new_m;
-   mesh_coarse = new_mc;
    FindNodePos();
 
    sfes = new FiniteElementSpace(mesh, new_fes->FEColl(), 1,
diff --git a/lib/vsvector3d.hpp b/lib/vsvector3d.hpp
index e8b1d75e..0efcffd6 100644
--- a/lib/vsvector3d.hpp
+++ b/lib/vsvector3d.hpp
@@ -33,6 +33,8 @@ class VisualizationSceneVector3d : public VisualizationSceneSolution3d
    Array<int> vflevel;
    Array<double> dvflevel;
 
+   int GetFunctionAutoRefineFactor() override;
+
 public:
    int ianim, ianimd, ianimmax, drawdisp;
 
@@ -44,12 +46,12 @@ class VisualizationSceneVector3d : public VisualizationSceneSolution3d
 
    virtual ~VisualizationSceneVector3d();
 
-   virtual std::string GetHelpString() const;
+   std::string GetHelpString() const override;
 
    void NPressed();
-   virtual void PrepareFlat();
-   virtual void Prepare();
-   virtual void PrepareLines();
+   void PrepareFlat() override;
+   void Prepare() override;
+   void PrepareLines() override;
 
    void PrepareFlat2();
    void PrepareLines2();
@@ -64,13 +66,13 @@ class VisualizationSceneVector3d : public VisualizationSceneSolution3d
    void SetScalarFunction();
    void ToggleScalarFunction();
 
-   virtual void PrepareCuttingPlane();
+   void PrepareCuttingPlane() override;
 
    void ToggleDisplacements() {drawdisp = (drawdisp+1)%2;};
 
-   virtual gl3::SceneInfo GetSceneObjs();
+   gl3::SceneInfo GetSceneObjs() override;
 
-   virtual void EventUpdateColors()
+   void EventUpdateColors() override
    { Prepare(); PrepareVectorField(); PrepareCuttingPlane(); };
 
    void ToggleVectorFieldLevel(int v);
diff --git a/tests/data b/tests/data
index 7457c82e..75f1a071 160000
--- a/tests/data
+++ b/tests/data
@@ -1 +1 @@
-Subproject commit 7457c82e272a521c7877ad54b1ac4e0533945377
+Subproject commit 75f1a0714b7176a3c9445bd83d2a90e52e1bde48