diff --git a/CHANGELOG b/CHANGELOG
index 60cedb5f..7d6fab7a 100644
--- a/CHANGELOG
+++ b/CHANGELOG
@@ -16,6 +16,10 @@ Version 4.2.1 (development)
 
 - Add support to visualize solutions on 1D elements embedded in 2D and 3D.
 
+- Added support for scalar integral finite elements, where calculation of the
+  surface normals was not implemented and was crashing GLVis. The normals are
+  approximately calculated from the point-wise projected value-based elements.
+
 - Added two new modes for visualization of vector fields in 2D, placing the
   arrows above the plotted surface and using a single color.
 
diff --git a/lib/vssolution.cpp b/lib/vssolution.cpp
index 0fc416a4..91e853eb 100644
--- a/lib/vssolution.cpp
+++ b/lib/vssolution.cpp
@@ -727,7 +727,50 @@ int VisualizationSceneSolution::GetRefinedValuesAndNormals(
       // In 1D we do not have well-defined normals.
       if (dim > 1)
       {
-         rsol->GetGradients(i, ir, tr);
+         const int map_type = rsol->FESpace()->GetFE(i)->GetMapType();
+         if (map_type == FiniteElement::MapType::VALUE)
+         {
+            rsol->GetGradients(i, ir, tr);
+         }
+         else if (map_type == FiniteElement::MapType::INTEGRAL)
+         {
+            FiniteElementSpace *fes = rsol->FESpace();
+            const FiniteElement *fe = fes->GetFE(i);
+            const int ndof = fe->GetDof();
+            const int ndim = fe->GetDim();
+            ElementTransformation *Trans = fes->GetElementTransformation(i);
+            DenseMatrix dshape(ndof, ndim);
+            Vector lval, gh(ndim), gcol;
+
+            rsol->GetElementDofValues(i, lval);
+
+            // Local projection to value-based FE
+            const IntegrationRule &nodes = fe->GetNodes();
+            for (int n = 0; n < nodes.GetNPoints(); n++)
+            {
+               const IntegrationPoint &ip = nodes.IntPoint(n);
+               Trans->SetIntPoint(&ip);
+               lval(n) /= Trans->Weight();//value = dof / |J|
+            }
+
+            // Gradient calculation
+            tr.SetSize(fe->GetDim(), ir.GetNPoints());
+            for (int q = 0; q < ir.GetNPoints(); q++)
+            {
+               const IntegrationPoint &ip = ir.IntPoint(q);
+               fe->CalcDShape(ip, dshape);
+               dshape.MultTranspose(lval, gh);
+               Trans->SetIntPoint(&ip);
+               tr.GetColumnReference(q, gcol);
+               const DenseMatrix &Jinv = Trans->InverseJacobian();
+               Jinv.MultTranspose(gh, gcol);
+            }
+         }
+         else
+         {
+            MFEM_ABORT("Unknown mapping type");
+         }
+
          normals.SetSize(3, tr.Width());
          for (int j = 0; j < tr.Width(); j++)
          {
diff --git a/lib/vssolution3d.cpp b/lib/vssolution3d.cpp
index 47c543ec..2a0daf52 100644
--- a/lib/vssolution3d.cpp
+++ b/lib/vssolution3d.cpp
@@ -3915,11 +3915,7 @@ void VisualizationSceneSolution3d::PrepareLevelSurf()
    {
       RefinedGeometry *RefG;
 #define GLVIS_SMOOTH_LEVELSURF_NORMALS
-#ifdef GLVIS_SMOOTH_LEVELSURF_NORMALS
-      const DenseMatrix *gp = &grad;
-#else
       const DenseMatrix *gp = NULL;
-#endif
 
       for (int ie = 0; ie < mesh->GetNE(); ie++)
       {
@@ -3928,7 +3924,52 @@ void VisualizationSceneSolution3d::PrepareLevelSurf()
          RefG = GLVisGeometryRefiner.Refine(geom, TimesToRefine);
          GridF->GetValues(ie, RefG->RefPts, vals, pointmat);
 #ifdef GLVIS_SMOOTH_LEVELSURF_NORMALS
-         GridF->GetGradients(ie, RefG->RefPts, grad);
+         const int map_type = GridF->FESpace()->GetFE(ie)->GetMapType();
+         if (map_type == FiniteElement::MapType::VALUE)
+         {
+            GridF->GetGradients(ie, RefG->RefPts, grad);
+            gp = &grad;
+         }
+         else if (map_type == FiniteElement::MapType::INTEGRAL)
+         {
+            FiniteElementSpace *fes = GridF->FESpace();
+            const FiniteElement *fe = fes->GetFE(ie);
+            const int ndof = fe->GetDof();
+            const int ndim = fe->GetDim();
+            ElementTransformation *Trans = fes->GetElementTransformation(ie);
+            const IntegrationRule &ir = RefG->RefPts;
+            DenseMatrix dshape(ndof, ndim);
+            Vector lval, gh(ndim), gcol;
+
+            GridF->GetElementDofValues(ie, lval);
+
+            // Local projection to value-based FE
+            const IntegrationRule &nodes = fe->GetNodes();
+            for (int n = 0; n < nodes.GetNPoints(); n++)
+            {
+               const IntegrationPoint &ip = nodes.IntPoint(n);
+               Trans->SetIntPoint(&ip);
+               lval(n) /= Trans->Weight();//value = dof / |J|
+            }
+
+            // Gradient calculation
+            grad.SetSize(fe->GetDim(), ir.GetNPoints());
+            for (int q = 0; q < ir.GetNPoints(); q++)
+            {
+               const IntegrationPoint &ip = ir.IntPoint(q);
+               fe->CalcDShape(ip, dshape);
+               dshape.MultTranspose(lval, gh);
+               Trans->SetIntPoint(&ip);
+               grad.GetColumnReference(q, gcol);
+               const DenseMatrix &Jinv = Trans->InverseJacobian();
+               Jinv.MultTranspose(gh, gcol);
+            }
+            gp = &grad;
+         }
+         else
+         {
+            MFEM_ABORT("Unknown mapping type");
+         }
 #endif
 
          Array<int> &RG = RefG->RefGeoms;
diff --git a/tests/data b/tests/data
index 81de802a..53fbee22 160000
--- a/tests/data
+++ b/tests/data
@@ -1 +1 @@
-Subproject commit 81de802ad63ba41cd1fcd3825fbb3d2471538ace
+Subproject commit 53fbee222cceb0581122cb72518bfc1b85766b4e