Merge pull request #285 from GLVis/higher-auto-ref

Higher auto-refinement limits

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.

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.

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;

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;

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()

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; }