Support for integral finite elements

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 FEs.
+
 
 Version 4.2 released on May 23, 2022
 ====================================

lib/vssolution.cpp

@@ -727,7 +727,44 @@ int VisualizationSceneSolution::GetRefinedValuesAndNormals(
       // In 1D we do not have well-defined normals.
       if (dim > 1)
       {
-         rsol->GetGradients(i, ir, tr);
+         if (rsol->FESpace()->GetFE(i)->GetMapType() == FiniteElement::VALUE)
+         {
+            rsol->GetGradients(i, ir, tr);
+         }
+         else
+         {
+            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();
+            }
+
+            // 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);
+            }
+         }
          normals.SetSize(3, tr.Width());
          for (int j = 0; j < tr.Width(); j++)
          {

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,47 @@ 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);
+         if (GridF->FESpace()->GetFE(ie)->GetMapType() == FiniteElement::VALUE)
+         {
+            GridF->GetGradients(ie, RefG->RefPts, grad);
+            gp = &grad;
+         }
+         else
+         {
+            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();
+            }
+
+            // 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;
+         }
 #endif
 
          Array<int> &RG = RefG->RefGeoms;