Fix test error

This reverts commit 9dbc125598.

doc/CodeDocumentation.conf.in

@@ -980,6 +980,7 @@ INPUT                  = @MFEM_SOURCE_DIR@/doc/CodeDocumentation.dox \
                          @MFEM_SOURCE_DIR@/miniapps/mtop \
                          @MFEM_SOURCE_DIR@/miniapps/multidomain \
                          @MFEM_SOURCE_DIR@/miniapps/navier \
+                         @MFEM_SOURCE_DIR@/miniapps/stabilized \
                          @MFEM_SOURCE_DIR@/miniapps/nurbs \
                          @MFEM_SOURCE_DIR@/miniapps/parelag \
                          @MFEM_SOURCE_DIR@/miniapps/performance \

examples/ex23.cpp

@@ -44,7 +44,7 @@ protected:
    BilinearForm *M;
    BilinearForm *K;
 
-   SparseMatrix Mmat, Kmat, Kmat0;
+   SparseMatrix Mmat, Kmat;
    SparseMatrix *T; // T = M + dt K
    real_t current_dt;
 
@@ -83,25 +83,24 @@ WaveOperator::WaveOperator(FiniteElementSpace &f,
    : SecondOrderTimeDependentOperator(f.GetTrueVSize(), (real_t) 0.0),
      fespace(f), M(NULL), K(NULL), T(NULL), current_dt(0.0), z(height)
 {
-   const real_t rel_tol = 1e-8;
-
-   fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
-
+   // Assemble Laplace matrix
    c2 = new ConstantCoefficient(speed*speed);
-
    K = new BilinearForm(&fespace);
    K->AddDomainIntegrator(new DiffusionIntegrator(*c2));
    K->Assemble();
 
-   Array<int> dummy;
-   K->FormSystemMatrix(dummy, Kmat0);
-   K->FormSystemMatrix(ess_tdof_list, Kmat);
-
+   // Assemble Mass matrix
    M = new BilinearForm(&fespace);
    M->AddDomainIntegrator(new MassIntegrator());
    M->Assemble();
+
+   // Apply Bcs
+   fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   K->FormSystemMatrix(ess_tdof_list, Kmat);
    M->FormSystemMatrix(ess_tdof_list, Mmat);
 
+   // Configure preconditioner
+   const real_t rel_tol = 1e-8;
    M_solver.iterative_mode = false;
    M_solver.SetRelTol(rel_tol);
    M_solver.SetAbsTol(0.0);
@@ -110,14 +109,13 @@ WaveOperator::WaveOperator(FiniteElementSpace &f,
    M_solver.SetPreconditioner(M_prec);
    M_solver.SetOperator(Mmat);
 
+   // Configure solver
    T_solver.iterative_mode = false;
    T_solver.SetRelTol(rel_tol);
    T_solver.SetAbsTol(0.0);
    T_solver.SetMaxIter(100);
    T_solver.SetPrintLevel(0);
    T_solver.SetPreconditioner(T_prec);
-
-   T = NULL;
 }
 
 void WaveOperator::Mult(const Vector &u, const Vector &du_dt,
@@ -126,9 +124,11 @@ void WaveOperator::Mult(const Vector &u, const Vector &du_dt,
    // Compute:
    //    d2udt2 = M^{-1}*-K(u)
    // for d2udt2
-   Kmat.Mult(u, z);
+   K->FullMult(u, z);
    z.Neg(); // z = -z
+   z.SetSubVector(ess_tdof_list, 0.0);
    M_solver.Mult(z, d2udt2);
+   d2udt2.SetSubVector(ess_tdof_list, 0.0);
 }
 
 void WaveOperator::ImplicitSolve(const real_t fac0, const real_t fac1,
@@ -142,14 +142,11 @@ void WaveOperator::ImplicitSolve(const real_t fac0, const real_t fac1,
       T = Add(1.0, Mmat, fac0, Kmat);
       T_solver.SetOperator(*T);
    }
-   Kmat0.Mult(u, z);
+   K->FullMult(u, z);
    z.Neg();
-
-   for (int i = 0; i < ess_tdof_list.Size(); i++)
-   {
-      z[ess_tdof_list[i]] = 0.0;
-   }
+   z.SetSubVector(ess_tdof_list, 0.0);
    T_solver.Mult(z, d2udt2);
+   d2udt2.SetSubVector(ess_tdof_list, 0.0);
 }
 
 void WaveOperator::SetParameters(const Vector &u)
@@ -314,7 +311,6 @@ int main(int argc, char *argv[])
          ess_bdr = 0;
       }
    }
-
    WaveOperator oper(fespace, ess_bdr, speed);
 
    u_gf.SetFromTrueDofs(u);

fem/bilininteg.cpp

@@ -1539,6 +1539,351 @@ const IntegrationRule &ConvectionIntegrator::GetRule(
    return GetRule(el,el,Trans);
 }
 
+
+void LaplaceIntegrator::AssembleElementMatrix(const FiniteElement &el,
+                                              ElementTransformation &Trans,
+                                              DenseMatrix &elmat )
+{
+   int nd = el.GetDof();
+   real_t w;
+
+   elmat.SetSize(nd);
+   shape.SetSize(nd);
+   laplace.SetSize(nd);
+
+   const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, el, Trans);
+
+   elmat = 0.0;
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      Trans.SetIntPoint (&ip);
+
+      el.CalcPhysShape(Trans, shape);
+      el.CalcPhysLaplacian(Trans, laplace);
+
+      w = Trans.Weight() * ip.weight * alpha;
+      if (Q)
+      {
+         w *= Q -> Eval(Trans, ip);
+      }
+      shape *= w;
+      AddMultVWt(shape, laplace, elmat);
+   }
+}
+
+void LaplaceIntegrator::AssembleElementMatrix2(const FiniteElement &trial_fe,
+                                               const FiniteElement &test_fe,
+                                               ElementTransformation &Trans,
+                                               DenseMatrix &elmat)
+{
+   int tr_nd = trial_fe.GetDof();
+   int te_nd = test_fe.GetDof();
+
+   elmat.SetSize(te_nd, tr_nd);
+   laplace.SetSize(tr_nd);
+   shape.SetSize(te_nd);
+
+   const IntegrationRule *ir = IntRule ? IntRule : &GetRule(trial_fe, test_fe,
+                                                            Trans);
+
+   elmat = 0.0;
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      Trans.SetIntPoint (&ip);
+
+      test_fe.CalcPhysShape(Trans, shape);
+      trial_fe.CalcPhysLaplacian(Trans, laplace);
+
+      real_t w = Trans.Weight() * ip.weight * alpha;
+      if (Q)
+      {
+         w *= Q -> Eval(Trans, ip);
+      }
+      AddMult_a_VWt(w, shape, laplace, elmat);
+   }
+}
+
+const IntegrationRule &LaplaceIntegrator::GetRule(
+   const FiniteElement &trial_fe,
+   const FiniteElement &test_fe,
+   ElementTransformation &Trans)
+{
+   int order = trial_fe.GetOrder() + test_fe.GetOrder();
+   return IntRules.Get(trial_fe.GetGeomType(), order);
+}
+
+void LaplaceGradIntegrator::AssembleElementMatrix(const FiniteElement &el,
+                                                  ElementTransformation &Trans,
+                                                  DenseMatrix &elmat )
+{
+   int nd = el.GetDof();
+   dim = el.GetDim();
+
+   elmat.SetSize(nd);
+   dshape.SetSize(nd,dim);
+   adjJ.SetSize(dim);
+   laplace.SetSize(nd);
+   vec2.SetSize(dim);
+   BdFidxT.SetSize(nd);
+
+   Vector vec1;
+
+   const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, el, Trans);
+
+   Q->Eval(Q_ir, Trans, *ir);
+
+   elmat = 0.0;
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      el.CalcDShape(ip, dshape);
+      el.CalcPhysLaplacian(Trans, laplace);
+
+      Trans.SetIntPoint(&ip);
+      CalcAdjugate(Trans.Jacobian(), adjJ);
+      Q_ir.GetColumnReference(i, vec1);
+      vec1 *= alpha * ip.weight;
+
+      adjJ.Mult(vec1, vec2);
+      dshape.Mult(vec2, BdFidxT);
+
+      AddMultVWt(BdFidxT, laplace, elmat);
+   }
+}
+
+void LaplaceGradIntegrator::AssembleElementMatrix2(const FiniteElement
+                                                   &trial_fe,
+                                                   const FiniteElement &test_fe,
+                                                   ElementTransformation &Trans,
+                                                   DenseMatrix &elmat)
+{
+   dim = trial_fe.GetDim();
+   int tr_nd = trial_fe.GetDof();
+   int te_nd = test_fe.GetDof();
+
+   elmat.SetSize(te_nd, tr_nd);
+   laplace.SetSize(tr_nd);
+   dshape.SetSize(te_nd,dim);
+   adjJ.SetSize(dim);
+   vec2.SetSize(dim);
+   BdFidxT.SetSize(te_nd);
+
+   Vector vec1;
+
+   const IntegrationRule *ir = IntRule ? IntRule : &GetRule(trial_fe, test_fe,
+                                                            Trans);
+
+   Q->Eval(Q_ir, Trans, *ir);
+
+   elmat = 0.0;
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      test_fe.CalcDShape(ip, dshape);
+      trial_fe.CalcPhysLaplacian(Trans, laplace);
+
+      Trans.SetIntPoint(&ip);
+      CalcAdjugate(Trans.Jacobian(), adjJ);
+      Q_ir.GetColumnReference(i, vec1);
+      vec1 *= alpha * ip.weight;
+
+      adjJ.Mult(vec1, vec2);
+      dshape.Mult(vec2, BdFidxT);
+
+      AddMultVWt(BdFidxT, laplace,elmat);
+   }
+}
+
+const IntegrationRule &LaplaceGradIntegrator::GetRule(
+   const FiniteElement &trial_fe,
+   const FiniteElement &test_fe,
+   ElementTransformation &Trans)
+{
+   int order = trial_fe.GetOrder() + test_fe.GetOrder();
+   return IntRules.Get(trial_fe.GetGeomType(), order);
+}
+
+void LaplaceLaplaceIntegrator::AssembleElementMatrix(const FiniteElement &el,
+                                                     ElementTransformation &Trans,
+                                                     DenseMatrix &elmat )
+{
+   int nd = el.GetDof();
+   real_t w;
+
+   elmat.SetSize(nd);
+   laplace.SetSize(nd);
+
+   const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, el, Trans);
+
+   elmat = 0.0;
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      Trans.SetIntPoint (&ip);
+
+      el.CalcPhysLaplacian(Trans, laplace);
+
+      w = Trans.Weight() * ip.weight * alpha;
+      if (Q)
+      {
+         w *= Q -> Eval(Trans, ip);
+      }
+      AddMult_a_VVt(w, laplace, elmat);
+   }
+}
+
+void LaplaceLaplaceIntegrator::AssembleElementMatrix2(const FiniteElement
+                                                      &trial_fe,
+                                                      const FiniteElement &test_fe,
+                                                      ElementTransformation &Trans,
+                                                      DenseMatrix &elmat)
+{
+   int
+   dim = trial_fe.GetDim();
+   int tr_nd = trial_fe.GetDof();
+   int te_nd = test_fe.GetDof();
+   real_t w;
+
+   elmat.SetSize(te_nd, tr_nd);
+   laplace.SetSize(tr_nd);
+   te_laplace.SetSize(te_nd);
+
+   const IntegrationRule *ir = IntRule ? IntRule : &GetRule(trial_fe, test_fe,
+                                                            Trans);
+
+   elmat = 0.0;
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      Trans.SetIntPoint (&ip);
+
+      trial_fe.CalcPhysLaplacian(Trans, laplace);
+      test_fe.CalcPhysLaplacian(Trans, te_laplace);
+
+      w = Trans.Weight() * ip.weight * alpha;
+      if (Q)
+      {
+         w *= Q -> Eval(Trans, ip);
+      }
+      AddMult_a_VWt(w, te_laplace, laplace, elmat);
+   }
+}
+
+const IntegrationRule &LaplaceLaplaceIntegrator::GetRule(
+   const FiniteElement &trial_fe,
+   const FiniteElement &test_fe,
+   ElementTransformation &Trans)
+{
+   int order = trial_fe.GetOrder() + test_fe.GetOrder();
+   return IntRules.Get(trial_fe.GetGeomType(), order);
+}
+
+void InverseEstimateIntegrator::AssembleElementMatrix(const FiniteElement &el,
+                                                      ElementTransformation &Trans,
+                                                      DenseMatrix &elmat )
+{
+   elmat = 0.0;
+
+   int nd = el.GetDof();
+   int dim = el.GetDim();
+
+   shape.SetSize(nd);
+   dshape.SetSize(nd,dim);
+   laplace.SetSize(nd);
+
+   lapmat.SetSize(nd,nd);
+   bimat.SetSize(nd,nd);
+   ovec.SetSize(nd);
+
+   real_t w,q;
+
+   const IntegrationRule *ir = IntRule;
+   if (ir == NULL)
+   {
+      int order = Trans.OrderGrad(&el) + Trans.Order() + el.GetOrder();
+      ir = &IntRules.Get(el.GetGeomType(), order);
+   }
+
+   bimat = 0.0;
+   lapmat = 0.0;
+   ovec = 0.0;
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      Trans.SetIntPoint(&ip);
+      w = Trans.Weight()*ip.weight;
+      if (Q)
+      {
+         q = Q->Eval(Trans, ip);
+      }
+
+      el.CalcPhysDShape(Trans, dshape);
+      AddMult_a_AAt(w*q, dshape, lapmat);
+
+      el.CalcPhysLaplacian(Trans, laplace);
+      AddMult_a_VVt(w*q*q, laplace, bimat);
+
+      el.CalcPhysShape(Trans, shape);
+      ovec.Add(w, shape);
+   }
+
+   // Power method
+   Vector x(nd);
+   x.Randomize(696383532);
+
+   // Correct nullspace + inverse
+   AddMult_a_VVt(1.0, ovec, lapmat);
+   DenseMatrixInverse L_inv(lapmat);
+
+   // DenseMatrix M_i, Q_i;
+   real_t alpha= 0.0, eval_i = 0.0, eval_prev = 0.0;
+
+   // Inverse power method
+   Vector x_tmp(nd);
+   int iter = 0;
+   const real_t rel_tol = 1e-4;
+
+   alpha = ovec*ovec;
+   ovec *= 1.0/sqrt(alpha);
+   do
+   {
+      // Othogonalize
+      alpha = x*ovec;
+      x.Add(-alpha, ovec);
+
+      // MatVec (2x)
+      bimat.Mult(x, x_tmp);
+      L_inv.Mult(x_tmp, x);
+
+      eval_prev = eval_i;
+      eval_i = x.Norml2();
+      x *= 1.0/eval_i;
+      ++iter;
+   }
+   while ((iter < 10000) && (fabs(eval_i - eval_prev)/fabs(eval_i) > rel_tol));
+   MFEM_VERIFY(fabs(eval_i - eval_prev)/fabs(eval_i) <= rel_tol,
+               "Inverse power method did not converge."
+               << "\n\t iter      = " << iter
+               << "\n\t eval_i    = " << eval_i
+               << "\n\t eval_prev = " << eval_prev
+               << "\n\t fabs(eval_i - eval_prev)/fabs(eval_i) = "
+               << fabs(eval_i - eval_prev)/fabs(eval_i));
+   cout<<"evev = "<<eval_i<<" "<<iter<<endl;
+}
+
+const IntegrationRule &InverseEstimateIntegrator::GetRule(
+   const FiniteElement &trial_fe,
+   const FiniteElement &test_fe,
+   ElementTransformation &Trans)
+{
+   // int order = Trans.OrderGrad(&trial_fe) + Trans.Order() + test_fe.GetOrder() - 2;
+   int order = trial_fe.GetOrder() + test_fe.GetOrder();
+   return IntRules.Get(trial_fe.GetGeomType(), order);
+}
+
+
 void VectorMassIntegrator::AssembleElementMatrix
 ( const FiniteElement &el, ElementTransformation &Trans,
   DenseMatrix &elmat )

fem/bilininteg.hpp

@@ -2450,6 +2450,135 @@ public:
                                       DenseMatrix &);
 };
 
+
+/// $\alpha (Q \Delta u, v)$
+class LaplaceIntegrator : public BilinearFormIntegrator
+{
+protected:
+   Coefficient *Q;
+   real_t alpha;
+
+private:
+   Vector laplace, shape;
+
+public:
+   LaplaceIntegrator(Coefficient &q, real_t a = 1.0)
+      : Q(&q) { alpha = a; }
+
+   virtual void AssembleElementMatrix(const FiniteElement &,
+                                      ElementTransformation &,
+                                      DenseMatrix &);
+
+   virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
+                                       const FiniteElement &test_fe,
+                                       ElementTransformation &Trans,
+                                       DenseMatrix &elmat);
+
+   static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
+                                         const FiniteElement &test_fe,
+                                         ElementTransformation &Trans);
+};
+
+/// $\alpha (u, Q \Delta v)$
+class TransposeLaplaceIntegrator : public TransposeIntegrator
+{
+public:
+   TransposeLaplaceIntegrator (Coefficient &q, real_t a = 1.0)
+      : TransposeIntegrator(new LaplaceIntegrator(q, a)) { }
+};
+
+/// $\alpha (\Delta u, Q \cdot \nabla v)$
+class LaplaceGradIntegrator : public BilinearFormIntegrator
+{
+protected:
+   VectorCoefficient *Q;
+   real_t alpha;
+   int dim;
+
+private:
+   Vector laplace, vec2, BdFidxT;
+   DenseMatrix dshape, adjJ, Q_ir;
+
+public:
+   LaplaceGradIntegrator(VectorCoefficient &q, real_t a = 1.0)
+      : Q(&q) { alpha = a; }
+
+   virtual void AssembleElementMatrix(const FiniteElement &,
+                                      ElementTransformation &,
+                                      DenseMatrix &);
+
+   virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
+                                       const FiniteElement &test_fe,
+                                       ElementTransformation &Trans,
+                                       DenseMatrix &elmat);
+
+   static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
+                                         const FiniteElement &test_fe,
+                                         ElementTransformation &Trans);
+};
+
+/// $\alpha (Q \cdot \nabla u,  \Delta v)$
+class GradLaplaceIntegrator : public TransposeIntegrator
+{
+public:
+   GradLaplaceIntegrator(VectorCoefficient &q, real_t a = 1.0)
+      : TransposeIntegrator(new LaplaceGradIntegrator(q, a)) { }
+};
+
+/// $\alpha (Q \Delta u, \Delta v)$
+class  LaplaceLaplaceIntegrator : public BilinearFormIntegrator
+{
+protected:
+   Coefficient *Q;
+   real_t alpha;
+
+private:
+   Vector laplace, te_laplace;
+
+public:
+   LaplaceLaplaceIntegrator(Coefficient &q, real_t a = 1.0)
+      : Q(&q) { alpha = a; }
+
+   virtual void AssembleElementMatrix(const FiniteElement &,
+                                      ElementTransformation &,
+                                      DenseMatrix &);
+
+   virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
+                                       const FiniteElement &test_fe,
+                                       ElementTransformation &Trans,
+                                       DenseMatrix &elmat);
+
+   static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
+                                         const FiniteElement &test_fe,
+                                         ElementTransformation &Trans);
+};
+
+// Alias for @LaplaceLaplaceIntegrator.
+using BiHarmonicIntegrator = LaplaceLaplaceIntegrator;
+
+/// Get the inverse estimate
+class InverseEstimateIntegrator : public BilinearFormIntegrator
+{
+protected:
+   Coefficient *Q;
+
+private:
+   Vector laplace, shape, ovec;//, vec2, BdFidxT;
+   DenseMatrix dshape, lapmat, bimat;//, adjJ, Q_ir;
+
+public:
+   InverseEstimateIntegrator(Coefficient &q)
+      : Q(&q) { }
+
+   virtual void AssembleElementMatrix(const FiniteElement &,
+                                      ElementTransformation &,
+                                      DenseMatrix &);
+
+   static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
+                                         const FiniteElement &test_fe,
+                                         ElementTransformation &Trans);
+};
+
 /** Class for integrating the bilinear form $a(u,v) := (Q u, v)$,
     where $u=(u_1,\dots,u_n)$ and $v=(v_1,\dots,v_n)$, $u_i$ and $v_i$ are defined
     by scalar FE through standard transformation. */

fem/coefficient.cpp

@@ -1504,6 +1504,295 @@ void CrossCrossCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
    M *= ((a == NULL ) ? aConst : a->Eval(T, ip) );
 }
 
+
+InverseEstimateCoefficient::InverseEstimateCoefficient(FiniteElementSpace *f)
+   : fes(f), Q(NULL), ir(NULL)
+{
+   ComputeInverseEstimates();
+}
+
+InverseEstimateCoefficient::InverseEstimateCoefficient(FiniteElementSpace *f,
+                                                       Coefficient &q)
+   : fes(f), Q(&q), ir(NULL)
+{
+   ComputeInverseEstimates();
+}
+
+GridFunction *InverseEstimateCoefficient::GetGridFunction()
+{
+   FiniteElementCollection* fec_ec = new L2_FECollection(0,
+                                                         fes ->GetMesh()->Dimension());
+   FiniteElementSpace *fes_ec = new FiniteElementSpace(fes ->GetMesh(), fec_ec);
+   GridFunction *gf = new GridFunction(fes_ec, elemInvEst.GetData());
+   gf->MakeOwner(fec_ec);
+   return gf;
+}
+
+void InverseEstimateCoefficient::ComputeInverseEstimates()
+{
+   elemInvEst.SetSize(fes -> GetNE());
+   SetIntRule(*fes->GetFE(0));
+   for (int i = 0; i < fes -> GetNE(); i++)
+   {
+      elemInvEst[i] = ElementInverseEstimate(*fes->GetFE(i),
+                                             *fes->GetElementTransformation(i));
+   }
+}
+
+void InverseEstimateCoefficient::SetIntRule(const FiniteElement &el)
+{
+   ir = &IntRules.Get(el.GetGeomType(), 2*el.GetOrder());
+}
+
+real_t InverseEstimateCoefficient::ElementInverseEstimate(
+   const FiniteElement &el,
+   ElementTransformation &Trans)
+{
+   if (el.GetOrder() < 2)
+   {
+      return std::numeric_limits<real_t>::min();
+   }
+
+   int nd = el.GetDof();
+   int dim = el.GetDim();
+
+   shape.SetSize(nd);
+   dshape.SetSize(nd,dim);
+   laplace.SetSize(nd);
+
+   lapmat.SetSize(nd,nd);
+   bimat.SetSize(nd,nd);
+   ovec.SetSize(nd);
+
+   real_t w,q = 1.0;
+
+   bimat = 0.0;
+   lapmat = 0.0;
+   ovec = 0.0;
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      Trans.SetIntPoint(&ip);
+      w = Trans.Weight()*ip.weight;
+      if (Q)
+      {
+         q = Q->Eval(Trans, ip);
+      }
+
+      el.CalcPhysDShape(Trans, dshape);
+      AddMult_a_AAt(w*q, dshape, lapmat);
+
+      el.CalcPhysLaplacian(Trans, laplace);
+      AddMult_a_VVt(w*q*q, laplace, bimat);
+
+      el.CalcPhysShape(Trans, shape);
+      ovec.Add(w, shape);
+   }
+   ovec *= 1.0/ovec.Norml2();
+
+   // Correct nullspace
+   AddMultVVt(ovec, lapmat);
+
+   // Return largest eigenvalue
+   return bimat.Eigenvalue(lapmat);
+}
+
+ElasticInverseEstimateCoefficient
+::ElasticInverseEstimateCoefficient(FiniteElementSpace *f)
+   : fes(f), Q(NULL), ir(NULL)
+{
+   ComputeInverseEstimates();
+}
+
+ElasticInverseEstimateCoefficient
+::ElasticInverseEstimateCoefficient(FiniteElementSpace *f,
+                                    Coefficient &q)
+   : fes(f), Q(&q), ir(NULL)
+{
+   ComputeInverseEstimates();
+}
+
+GridFunction *ElasticInverseEstimateCoefficient::GetGridFunction()
+{
+   FiniteElementCollection* fec_ec = new L2_FECollection(0,
+                                                         fes ->GetMesh()->Dimension());
+   FiniteElementSpace *fes_ec = new FiniteElementSpace(fes ->GetMesh(), fec_ec);
+   GridFunction *gf = new GridFunction(fes_ec, elemInvEst.GetData());
+   gf->MakeOwner(fec_ec);
+   return gf;
+}
+
+void ElasticInverseEstimateCoefficient::ComputeInverseEstimates()
+{
+   elemInvEst.SetSize(fes -> GetNE());
+   SetIntRule(*fes->GetFE(0));
+   int dim = fes->GetFE(0)->GetDim();
+
+   emat.SetSize(dim,dim);
+   divmat.SetSize(dim,dim);
+   for (int i = 0; i < dim; i++)
+   {
+      for (int j = 0; j < dim; j++)
+      {
+         emat(i,j)= new DenseMatrix();
+         divmat(i,j)= new DenseMatrix();
+      }
+   }
+
+   hmap.SetSize(dim,dim);
+
+   if (dim == 2)
+   {
+      hmap(0,0) = 0;
+      hmap(0,1) =  hmap(1,0) =  1;
+      hmap(1,1) = 2;
+   }
+   else if (dim == 2)
+   {
+      hmap(0,0) = 0;
+      hmap(0,1) = hmap(1,0) = 1;
+      hmap(0,2) = hmap(2,0) = 2;
+      hmap(1,1) = 3;
+      hmap(1,2) = hmap(2,1) = 4;
+      hmap(2,2) = 5;
+   }
+   else
+   {
+      mfem_error("Only implemented for 2D and 3D");
+   }
+
+   for (int i = 0; i < fes -> GetNE(); i++)
+   {
+      elemInvEst[i] = ElementInverseEstimate(*fes->GetFE(i),
+                                             *fes->GetElementTransformation(i));
+   }
+}
+
+void ElasticInverseEstimateCoefficient::SetIntRule(const FiniteElement &el)
+{
+   ir = &IntRules.Get(el.GetGeomType(), 2*el.GetOrder());
+}
+
+real_t ElasticInverseEstimateCoefficient::ElementInverseEstimate(
+   const FiniteElement &el,
+   ElementTransformation &Trans)
+{
+   // if (el.GetDerivType() != (int) FiniteElement::HESS)
+   //  {
+   //     return std::numeric_limits<real_t>::min();
+   //  }
+
+   int nd = el.GetDof();
+   int dim = el.GetDim();
+
+   shape.SetSize(nd);
+   dshape.SetSize(nd,dim);
+   hshape.SetSize(nd,dim*(dim+1)/2);
+
+   for (int i = 0; i < dim; i++)
+   {
+      for (int j = 0; j < dim; j++)
+      {
+         emat(i,j)->SetSize(nd,nd);
+         *emat(i,j) = 0.0;
+         divmat(i,j)->SetSize(nd,nd);
+         *divmat(i,j) = 0.0;
+      }
+   }
+
+   real_t w,q = 1.0;
+   for (int ii = 0; ii < ir->GetNPoints(); ii++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(ii);
+      Trans.SetIntPoint(&ip);
+      w = Trans.Weight()*ip.weight;
+      if (Q)
+      {
+         q = Q->Eval(Trans, ip);
+      }
+
+      el.CalcPhysDShape(Trans, dshape);
+      for (int i = 0; i < dim; i++)
+      {
+         for (int j = 0; j < dim; j++)
+         {
+            AddMult_a_VVt(w*q, Vector(dshape.GetColumn(i),nd), *emat(j,j));
+
+            AddMult_a_VWt(w*q, Vector(dshape.GetColumn(i),nd),
+                          Vector(dshape.GetColumn(j),nd), *emat(j,i));
+
+            AddMult_a_VWt(w*q, Vector(dshape.GetColumn(j),nd),
+                          Vector(dshape.GetColumn(i),nd), *emat(i,j));
+
+            AddMult_a_VVt(w*q, Vector(dshape.GetColumn(j),nd), *emat(i,i));
+         }
+      }
+
+      el.CalcPhysHessian(Trans, hshape);
+      for (int i = 0; i < dim; i++)
+      {
+         for (int j = 0; j < dim; j++)
+         {
+            for (int k = 0; k < dim; k++)
+            {
+               AddMult_a_VWt(w*q*q, Vector(hshape.GetColumn(hmap(i,i)),nd),
+                             Vector(hshape.GetColumn(hmap(k,k)),nd), *divmat(j,j));
+
+               AddMult_a_VWt(w*q*q, Vector(hshape.GetColumn(hmap(i,i)),nd),
+                             Vector(hshape.GetColumn(hmap(k,j)),nd), *divmat(j,k));
+
+               AddMult_a_VWt(w*q*q, Vector(hshape.GetColumn(hmap(i,j)),nd),
+                             Vector(hshape.GetColumn(hmap(k,k)),nd), *divmat(i,j));
+
+               AddMult_a_VWt(w*q*q, Vector(hshape.GetColumn(hmap(i,j)),nd),
+                             Vector(hshape.GetColumn(hmap(k,j)),nd), *divmat(i,k));
+            }
+         }
+      }
+   }
+
+   // Collect matrices
+   emat_tot.SetSize(nd*dim,nd*dim);
+   divmat_tot.SetSize(nd*dim,nd*dim);
+   for (int i = 0; i < dim; i++)
+   {
+      for (int j = 0; j < dim; j++)
+      {
+         emat_tot  .SetSubMatrix(i*nd, j*nd, *emat(i,j));
+         divmat_tot.SetSubMatrix(i*nd, j*nd, *divmat(i,j));
+      }
+   }
+
+   // Correct nullspace
+   DenseMatrix ns;
+   emat_tot.NullSpace(ns, 1e-10);
+   for (int i = 0; i < ns.Width(); i++)
+   {
+      AddMultVVt(Vector(ns.GetColumn(i),nd*dim), emat_tot);
+   }
+
+   // Return largest eigenvalue
+   return divmat_tot.Eigenvalue(emat_tot);
+}
+
+ElasticInverseEstimateCoefficient::~ElasticInverseEstimateCoefficient()
+{
+   for (int i = 0; i < emat.NumRows(); i++)
+   {
+      for (int j = 0; j < emat.NumCols(); j++)
+      {
+         delete emat(i,j);
+      }
+   }
+   for (int i = 0; i < divmat.NumRows(); i++)
+   {
+      for (int j = 0; j < divmat.NumCols(); j++)
+      {
+         delete divmat(i,j);
+      }
+   }
+}
+
 real_t LpNormLoop(real_t p, Coefficient &coeff, Mesh &mesh,
                   const IntegrationRule *irs[])
 {

fem/coefficient.hpp

@@ -2328,6 +2328,127 @@ public:
 };
 ///@}
 
+/** @brief
+*/
+class InverseEstimateCoefficient : public Coefficient
+{
+private:
+   ///
+   Vector elemInvEst;
+   /// FE space on which the grid function lives. Owned if #fec is not NULL.
+   FiniteElementSpace *fes;
+
+   ///
+   const IntegrationRule *ir;
+
+   ///
+   Coefficient *Q;
+   Vector laplace, shape, ovec, evec;
+   DenseMatrix dshape, lapmat, bimat;
+
+   ///
+   void SetIntRule(const FiniteElement &el);
+
+   ///
+   void ComputeInverseEstimates();
+
+   real_t ElementInverseEstimate(const FiniteElement &el,
+                                 ElementTransformation &Trans);
+
+public:
+   ///
+   InverseEstimateCoefficient(FiniteElementSpace *f);
+   InverseEstimateCoefficient(FiniteElementSpace *f, Coefficient &q);
+
+   /// Caller gets owner ship of GridFunction and
+   GridFunction *GetGridFunction();
+
+   /// Reset the scalar factor
+   void SetDiffusion(Coefficient &q)
+   {
+      if (Q != &q)
+      {
+         Q = &q;
+         ComputeInverseEstimates();
+      }
+   }
+   /// Return the scalar factor
+   Coefficient * GetDiffusion() const { return Q; }
+
+   /// Evaluate the coefficient at @a ip.
+   virtual real_t Eval(ElementTransformation &T,
+                       const IntegrationPoint &ip)
+   { return elemInvEst[T.ElementNo]; }
+
+};
+
+class ElasticInverseEstimateCoefficient : public Coefficient
+{
+private:
+   ///
+   Vector elemInvEst;
+   /// FE space on which the grid function lives. Owned if #fec is not NULL.
+   FiniteElementSpace *fes;
+
+   ///
+   const IntegrationRule *ir;
+
+   ///
+   Coefficient *Q;
+   Vector shape, ovec, evec;
+   DenseMatrix dshape, hshape, emat_tot, divmat_tot;
+   Array2D<DenseMatrix*> emat,divmat;
+   Array2D<int> hmap;
+
+   ///
+   void SetIntRule(const FiniteElement &el);
+
+   ///
+   void ComputeInverseEstimates();
+
+   ///
+   real_t ElementInverseEstimate(const FiniteElement &el,
+                                 ElementTransformation &Trans);
+
+public:
+   ///
+   ElasticInverseEstimateCoefficient(FiniteElementSpace *f);
+   ElasticInverseEstimateCoefficient(FiniteElementSpace *f, Coefficient &q);
+
+   /// Caller gets owner ship of GridFunction and
+   GridFunction *GetGridFunction();
+
+   /// Reset the scalar factor
+   void SetDiffusion(Coefficient &q)
+   {
+      if (Q != &q)
+      {
+         Q = &q;
+         ComputeInverseEstimates();
+      }
+   }
+   void SetShearModulus(Coefficient &q) { SetDiffusion(q);}
+
+   /// Return the scalar factor
+   Coefficient * GetDiffusion() const { return Q; }
+   Coefficient * GetModulus() const { return GetDiffusion(); }
+
+   /// Evaluate the coefficient at @a ip.
+   virtual real_t Eval(ElementTransformation &T,
+                       const IntegrationPoint &ip)
+   { return elemInvEst[T.ElementNo]; }
+
+   // Destructor
+   ~ElasticInverseEstimateCoefficient();
+
+};
+
+
+///@}
+
+
+
+
 /** @brief Vector quadrature function coefficient which requires that the
     quadrature rules used for this vector coefficient be the same as those that
     live within the supplied QuadratureFunction. */

fem/fe/fe_base.cpp

@@ -221,7 +221,7 @@ void FiniteElement::CalcPhysLaplacian(ElementTransformation &Trans,
    {
       for (int nd = 0; nd < dof; nd++)
       {
-         Laplacian[nd] = hess(nd,0) + hess(nd,4) + hess(nd,5);
+         Laplacian[nd] = hess(nd,0) + hess(nd,3) + hess(nd,5);
       }
    }
    else if (dim == 2)
@@ -259,10 +259,10 @@ void FiniteElement::CalcPhysLinLaplacian(ElementTransformation &Trans,
       scale[1] = 2*Gij(0,1);
       scale[2] = 2*Gij(0,2);
 
-      scale[3] = 2*Gij(1,2);
-      scale[4] =   Gij(2,2);
+      scale[3] =   Gij(1,1);
+      scale[4] = 2*Gij(1,2);
 
-      scale[5] =   Gij(1,1);
+      scale[5] =   Gij(2,2);
    }
    else if (dim == 2)
    {
@@ -299,12 +299,12 @@ void  FiniteElement::CalcPhysHessian(ElementTransformation &Trans,
       map[2] = 2;
 
       map[3] = 1;
-      map[4] = 5;
-      map[5] = 3;
+      map[4] = 3;
+      map[5] = 4;
 
       map[6] = 2;
-      map[7] = 3;
-      map[8] = 4;
+      map[7] = 4;
+      map[8] = 5;
    }
    else if (dim == 2)
    {

fem/fe/fe_base.hpp

@@ -299,7 +299,8 @@ public:
       NONE, ///< No derivatives implemented
       GRAD, ///< Implements CalcDShape methods
       DIV,  ///< Implements CalcDivShape methods
-      CURL  ///< Implements CalcCurlShape methods
+      CURL, ///< Implements CalcCurlShape methods
+      HESS  ///< Implements CalcHessian & CalcDShape methods
    };
 
    /** @brief Construct FiniteElement with given
@@ -356,7 +357,7 @@ public:
 
    /** @brief Returns the FiniteElement::DerivType of the element describing the
        spatial derivative method implemented, one of {NONE, GRAD,
-       DIV, CURL}. */
+       DIV, CURL, HESS}. */
    int GetDerivType() const { return deriv_type; }
 
    /** @brief Returns the FiniteElement::DerivType of the element describing how

fem/fe/fe_nurbs.cpp

@@ -349,10 +349,10 @@ void NURBS3DFiniteElement::CalcHessian (const IntegrationPoint &ip,
             d2sum[1] += ( hessian(o,1) = dsx*dsy*sz*weights(o) );
             d2sum[2] += ( hessian(o,2) = dsx*sy*dsz*weights(o) );
 
-            d2sum[3] += ( hessian(o,3) = sx*dsy*dsz*weights(o) );
+            d2sum[3] += ( hessian(o,3) = sx*d2sy*sz*weights(o) );
+            d2sum[4] += ( hessian(o,4) = sx*dsy*dsz*weights(o) );
 
-            d2sum[4] += ( hessian(o,4) = sx*sy*d2sz*weights(o) );
-            d2sum[5] += ( hessian(o,5) = sx*d2sy*sz*weights(o) );
+            d2sum[5] += ( hessian(o,5) = sx*sy*d2sz*weights(o) );
          }
       }
    }
@@ -387,17 +387,17 @@ void NURBS3DFiniteElement::CalcHessian (const IntegrationPoint &ip,
                      + u[o]*sum*(2*dsum[0]*dsum[2] - d2sum[2]);
 
       hessian(o,3) = hessian(o,3)*sum
-                     - du(o,1)*sum*dsum[2]
-                     - du(o,2)*sum*dsum[1]
-                     + u[o]*sum*(2*dsum[1]*dsum[2] - d2sum[3]);
+                     - 2*du(o,1)*sum*dsum[1]
+                     + u[o]*sum*(2*dsum[1]*dsum[1] - d2sum[3]);
 
       hessian(o,4) = hessian(o,4)*sum
-                     - 2*du(o,2)*sum*dsum[2]
-                     + u[o]*sum*(2*dsum[2]*dsum[2] - d2sum[4]);
+                     - du(o,1)*sum*dsum[2]
+                     - du(o,2)*sum*dsum[1]
+                     + u[o]*sum*(2*dsum[1]*dsum[2] - d2sum[4]);
 
       hessian(o,5) = hessian(o,5)*sum
-                     - 2*du(o,1)*sum*dsum[1]
-                     + u[o]*sum*(2*dsum[1]*dsum[1] - d2sum[5]);
+                     - 2*du(o,2)*sum*dsum[2]
+                     + u[o]*sum*(2*dsum[2]*dsum[2] - d2sum[5]);
    }
 }
 

fem/gridfunc.cpp

@@ -3051,6 +3051,57 @@ real_t GridFunction::ComputeDivError(
    return (error < 0.0) ? -sqrt(-error) : sqrt(error);
 }
 
+real_t GridFunction::ComputeLaplaceError(
+   Coefficient *exlap, const IntegrationRule *irs[]) const
+{
+   real_t error = 0.0, a;
+   const FiniteElement *fe;
+   ElementTransformation *Tr;
+   Array<int> dofs;
+   int intorder, fdof;
+   Vector laplace;
+
+   for (int i = 0; i < fes->GetNE(); i++)
+   {
+      laplace.SetSize(fdof);
+      fe = fes->GetFE(i);
+      Tr = fes->GetElementTransformation(i);
+      intorder = 2*fe->GetOrder() + 3;
+      const IntegrationRule *ir;
+      if (irs)
+      {
+         ir = irs[fe->GetGeomType()];
+      }
+      else
+      {
+         ir = &(IntRules.Get(fe->GetGeomType(), intorder));
+      }
+      fes->GetElementDofs(i, dofs);
+      fdof = fe->GetDof();
+      laplace.SetSize(fdof);
+      for (int j = 0; j < ir->GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(j);
+         Tr->SetIntPoint(&ip);
+         fe->CalcPhysLaplacian(*Tr, laplace);
+         a = 0;
+         for (int k = 0; k < fdof; k++)
+            if (dofs[k] >= 0)
+            {
+               a += (*this)(dofs[k]) * laplace(k);
+            }
+            else
+            {
+               a -= (*this)(-1-dofs[k]) * laplace(k);
+            }
+         a -= exlap->Eval(*Tr, ip);
+         error += ip.weight * Tr->Weight() * a * a;
+      }
+   }
+
+   return (error < 0.0) ? -sqrt(-error) : sqrt(error);
+}
+
 real_t GridFunction::ComputeDGFaceJumpError(Coefficient *exsol,
                                             Coefficient *ell_coeff,
                                             class JumpScaling jump_scaling,

fem/gridfunc.hpp

@@ -533,6 +533,10 @@ public:
    virtual real_t ComputeDivError(Coefficient *exdiv,
                                   const IntegrationRule *irs[] = NULL) const;
 
+   /// Returns ||lap u_ex - lap u_h||_L2 for H1 elements
+   virtual real_t ComputeLaplaceError(Coefficient *exlap,
+                                      const IntegrationRule *irs[] = NULL) const;
+
    /// Returns the Face Jumps error for L2 elements. The error can be weighted
    /// by a constant nu, by nu/h, or nu*p^2/h, depending on the value of
    /// @a jump_scaling.

fem/lininteg.cpp

@@ -123,6 +123,35 @@ void DomainLFGradIntegrator::AssembleDeltaElementVect(
    dshape.Mult(Qvec, elvect);
 }
 
+void DomainLFLaplaceIntegrator::AssembleRHSElementVect(const FiniteElement &el,
+                                                       ElementTransformation &Tr,
+                                                       Vector &elvect)
+{
+   int dof = el.GetDof();
+
+   laplace.SetSize(dof);       // vector of size dof
+   elvect.SetSize(dof);
+   elvect = 0.0;
+
+   const IntegrationRule *ir = NULL;//IntRule;
+   if (ir == NULL)
+   {
+      ir = &IntRules.Get(el.GetGeomType(), oa * el.GetOrder() + ob + 4);
+   }
+
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+
+      Tr.SetIntPoint (&ip);
+      real_t val = Tr.Weight() * Q.Eval(Tr, ip) * alpha;
+
+      el.CalcPhysLaplacian(Tr, laplace);
+
+      add(elvect, ip.weight * val, laplace, elvect);
+   }
+}
+
 void BoundaryLFIntegrator::AssembleRHSElementVect(
    const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
 {

fem/lininteg.hpp

@@ -174,6 +174,29 @@ public:
    using LinearFormIntegrator::AssembleRHSElementVect;
 };
 
+/// Class for domain integrator $ L(v) := (f, \Delta v) $
+class DomainLFLaplaceIntegrator : public LinearFormIntegrator
+{
+private:
+   Vector laplace;
+   Coefficient &Q;
+   real_t alpha;
+   int oa, ob;
+public:
+   /// Constructs the domain integrator $ (Q, \nabla v) $
+   DomainLFLaplaceIntegrator(Coefficient &QF, real_t alp = 1.0, int a = 2,
+                             int b = 0)
+      : Q(QF), oa(a), ob(b)   { alpha = alp; }
+
+   /** Given a particular Finite Element and a transformation (Tr)
+       computes the element right hand side element vector, elvect. */
+   virtual void AssembleRHSElementVect(const FiniteElement &el,
+                                       ElementTransformation &Tr,
+                                       Vector &elvect);
+
+   using LinearFormIntegrator::AssembleRHSElementVect;
+};
+
 
 /// Class for boundary integration $ L(v) := (g, v) $
 class BoundaryLFIntegrator : public LinearFormIntegrator

linalg/densemat.cpp

@@ -32,7 +32,6 @@
 #define copysign _copysign
 #endif
 
-
 namespace mfem
 {
 
@@ -907,169 +906,6 @@ void DenseMatrix::FNorm(real_t &scale_factor, real_t &scaled_fnorm2) const
    scaled_fnorm2 = fnorm2;
 }
 
-void dsyevr_Eigensystem(DenseMatrix &a, Vector &ev, DenseMatrix *evect)
-{
-#ifdef MFEM_USE_LAPACK
-   ev.SetSize(a.Width());
-
-   char      JOBZ     = 'N';
-   char      RANGE    = 'A';
-   char      UPLO     = 'U';
-   int       N        = a.Width();
-   real_t   *A        = new real_t[N*N];
-   int       LDA      = N;
-   real_t    VL       = 0.0;
-   real_t    VU       = 1.0;
-   int       IL       = 0;
-   int       IU       = 1;
-   real_t    ABSTOL   = 0.0;
-   int       M;
-   real_t   *W        = ev.GetData();
-   real_t   *Z        = NULL;
-   int       LDZ      = 1;
-   int      *ISUPPZ   = new int[2*N];
-   int       LWORK    = -1; // query optimal (double) workspace size
-   real_t    QWORK;
-   real_t   *WORK     = NULL;
-   int       LIWORK   = -1; // query optimal (int) workspace size
-   int       QIWORK;
-   int      *IWORK    = NULL;
-   int       INFO;
-
-   if (evect) // Compute eigenvectors too
-   {
-      evect->SetSize(N);
-
-      JOBZ     = 'V';
-      Z        = evect->Data();
-      LDZ      = N;
-   }
-
-   int hw = a.Height() * a.Width();
-   real_t *data = a.Data();
-
-   for (int i = 0; i < hw; i++)
-   {
-      A[i] = data[i];
-   }
-
-   MFEM_LAPACK_PREFIX(syevr_)(&JOBZ, &RANGE, &UPLO, &N, A, &LDA, &VL, &VU, &IL,
-                              &IU, &ABSTOL, &M, W, Z, &LDZ, ISUPPZ, &QWORK,
-                              &LWORK, &QIWORK, &LIWORK, &INFO);
-
-   LWORK  = (int) QWORK;
-   LIWORK = QIWORK;
-
-   WORK  = new real_t[LWORK];
-   IWORK = new int[LIWORK];
-
-   MFEM_LAPACK_PREFIX(syevr_)(&JOBZ, &RANGE, &UPLO, &N, A, &LDA, &VL, &VU, &IL,
-                              &IU, &ABSTOL, &M, W, Z, &LDZ, ISUPPZ, WORK,
-                              &LWORK, IWORK, &LIWORK, &INFO);
-
-   if (INFO != 0)
-   {
-      mfem::err << "dsyevr_Eigensystem(...): DSYEVR error code: "
-                << INFO << endl;
-      mfem_error();
-   }
-
-#ifdef MFEM_DEBUG
-   if (M < N)
-   {
-      mfem::err << "dsyevr_Eigensystem(...):\n"
-                << " DSYEVR did not find all eigenvalues "
-                << M << "/" << N << endl;
-      mfem_error();
-   }
-   if (CheckFinite(W, N) > 0)
-   {
-      mfem_error("dsyevr_Eigensystem(...): inf/nan values in W");
-   }
-   if (CheckFinite(Z, N*N) > 0)
-   {
-      mfem_error("dsyevr_Eigensystem(...): inf/nan values in Z");
-   }
-   VU = 0.0;
-   for (IL = 0; IL < N; IL++)
-      for (IU = 0; IU <= IL; IU++)
-      {
-         VL = 0.0;
-         for (M = 0; M < N; M++)
-         {
-            VL += Z[M+IL*N] * Z[M+IU*N];
-         }
-         if (IU < IL)
-         {
-            VL = fabs(VL);
-         }
-         else
-         {
-            VL = fabs(VL-1.0);
-         }
-         if (VL > VU)
-         {
-            VU = VL;
-         }
-         if (VU > 0.5)
-         {
-            mfem::err << "dsyevr_Eigensystem(...):"
-                      << " Z^t Z - I deviation = " << VU
-                      << "\n W[max] = " << W[N-1] << ", W[min] = "
-                      << W[0] << ", N = " << N << endl;
-            mfem_error();
-         }
-      }
-   if (VU > 1e-9)
-   {
-      mfem::err << "dsyevr_Eigensystem(...):"
-                << " Z^t Z - I deviation = " << VU
-                << "\n W[max] = " << W[N-1] << ", W[min] = "
-                << W[0] << ", N = " << N << endl;
-   }
-   if (VU > 1e-5)
-   {
-      mfem_error("dsyevr_Eigensystem(...): ERROR: ...");
-   }
-   VU = 0.0;
-   for (IL = 0; IL < N; IL++)
-      for (IU = 0; IU < N; IU++)
-      {
-         VL = 0.0;
-         for (M = 0; M < N; M++)
-         {
-            VL += Z[IL+M*N] * W[M] * Z[IU+M*N];
-         }
-         VL = fabs(VL-data[IL+N*IU]);
-         if (VL > VU)
-         {
-            VU = VL;
-         }
-      }
-   if (VU > 1e-9)
-   {
-      mfem::err << "dsyevr_Eigensystem(...):"
-                << " max matrix deviation = " << VU
-                << "\n W[max] = " << W[N-1] << ", W[min] = "
-                << W[0] << ", N = " << N << endl;
-   }
-   if (VU > 1e-5)
-   {
-      mfem_error("dsyevr_Eigensystem(...): ERROR: ...");
-   }
-#endif
-
-   delete [] IWORK;
-   delete [] WORK;
-   delete [] ISUPPZ;
-   delete [] A;
-#else
-   MFEM_CONTRACT_VAR(a);
-   MFEM_CONTRACT_VAR(ev);
-   MFEM_CONTRACT_VAR(evect);
-#endif
-}
-
 void dsyev_Eigensystem(DenseMatrix &a, Vector &ev, DenseMatrix *evect)
 {
 #ifdef MFEM_USE_LAPACK
@@ -1127,23 +963,222 @@ void dsyev_Eigensystem(DenseMatrix &a, Vector &ev, DenseMatrix *evect)
 #endif
 }
 
-void DenseMatrix::Eigensystem(Vector &ev, DenseMatrix *evect)
+void dsyevr_Eigensystem(DenseMatrix &a, Vector &ev, DenseMatrix *evect,
+                        char RANGE, real_t VL, real_t VU, int IL, int IU)
 {
 #ifdef MFEM_USE_LAPACK
+   ev.SetSize(a.Width());
 
-   // dsyevr_Eigensystem(*this, ev, evect);
+   char      JOBZ     = 'N';
+   char      UPLO     = 'U';
+   int       N        = a.Width();
+   real_t   *A        = new real_t[N*N];
+   int       LDA      = N;
+   real_t    ABSTOL   = 0.0;
+   int       M;
+   real_t   *W        = ev.GetData();
+   real_t   *Z        = NULL;
+   int       LDZ      = 1;
+   int      *ISUPPZ   = new int[2*N];
+   int       LWORK    = -1; // query optimal (double) workspace size
+   real_t    QWORK;
+   real_t   *WORK     = NULL;
+   int       LIWORK   = -1; // query optimal (int) workspace size
+   int       QIWORK;
+   int      *IWORK    = NULL;
+   int       INFO;
 
-   dsyev_Eigensystem(*this, ev, evect);
+   if (evect) // Compute eigenvectors too
+   {
+      evect->SetSize(N);
 
+      JOBZ     = 'V';
+      Z        = evect->Data();
+      LDZ      = N;
+   }
+
+
+   int hw = a.Height() * a.Width();
+   real_t *data = a.Data();
+
+   for (int i = 0; i < hw; i++)
+   {
+      A[i] = data[i];
+   }
+
+   MFEM_LAPACK_PREFIX(syevr_)(&JOBZ, &RANGE, &UPLO, &N, A, &LDA, &VL, &VU, &IL,
+                              &IU, &ABSTOL, &M, W, Z, &LDZ, ISUPPZ, &QWORK,
+                              &LWORK, &QIWORK, &LIWORK, &INFO);
+
+   LWORK  = (int) QWORK;
+   LIWORK = QIWORK;
+
+   WORK  = new real_t[LWORK];
+   IWORK = new int[LIWORK];
+
+   MFEM_LAPACK_PREFIX(syevr_)(&JOBZ, &RANGE, &UPLO, &N, A, &LDA, &VL, &VU, &IL,
+                              &IU, &ABSTOL, &M, W, Z, &LDZ, ISUPPZ, WORK,
+                              &LWORK, IWORK, &LIWORK, &INFO);
+
+   if (INFO != 0)
+   {
+      mfem::err << "dsyevr_Eigensystem(...): DSYEVR error code: "
+                << INFO << endl;
+      mfem_error();
+   }
+
+
+   if (evect) // Compute eigenvectors too
+   {
+      evect->SetSize(N,M);
+   }
+
+
+#ifdef MFEM_DEBUG
+   for (int i = 0; i < N; i++)
+   {
+      for (int j = i+1; j < N; j++)
+      {
+         if (fabs(data[i+N*j]-data[j+N*i]) > 1e-9)
+         {
+            mfem::err << "dsyevr_Eigensystem(...): matrix not symmetric\n"
+                      << " data["<<i<<"+N*"<<j<<"] = " << data[i+N*j]<<endl
+                      << " data["<<j<<"+N*"<<i<<"] = " << data[j+N*i]<<endl;
+            mfem_error();
+         }
+      }
+   }
+
+   if ((M < N) && (RANGE == 'A'))
+   {
+      mfem::err << "dsyevr_Eigensystem(...):\n"
+                << " DSYEVR did not find all eigenvalues "
+                << M << "/" << N << endl;
+      mfem_error();
+   }
+
+   // Check eigen vectors
+   if (evect)
+   {
+      if (CheckFinite(Z, N*N) > 0)
+      {
+         mfem_error("dsyevr_Eigensystem(...): inf/nan values in Z");
+      }
+
+      // Check orthogonality of eigenvectors
+      real_t u = 0.0;
+      for (int i = 0; i < M; i++)
+      {
+         for (int j = 0; j <= i; j++)
+         {
+            real_t l = 0.0;
+            for (int k = 0; k < N; k++)
+            {
+               l += Z[k+i*N] * Z[k+j*N];
+            }
+            if (j < i)
+            {
+               l = fabs(l);
+            }
+            else
+            {
+               l = fabs(l-1.0);
+            }
+            if (l > u)
+            {
+               u = l;
+            }
+            if (u > 0.5)
+            {
+               mfem::err << "dsyevr_Eigensystem(...):"
+                         << " Z^t Z - I deviation = " << u
+                         << "\n W[max] = " << W[N-1] << ", W[min] = "
+                         << W[0] << ", N = " << N << endl;
+               mfem_error();
+            }
+         }
+      }
+      if (u > 1e-9)
+      {
+         mfem::err << "dsyevr_Eigensystem(...):"
+                   << " Z^t Z - I deviation = " << u
+                   << "\n W[max] = " << W[N-1] << ", W[min] = "
+                   << W[0] << ", N = " << N << endl;
+      }
+      if (u > 1e-5)
+      {
+         mfem_error("dsyevr_Eigensystem(...): ERROR: ...");
+      }
+   }
+
+   // Check eigenvalues
+   if (CheckFinite(W, M) > 0)
+   {
+      mfem_error("dsyevr_Eigensystem(...): inf/nan values in W");
+   }
+
+   // Check if eigen decomposition generates original matrix
+   if (evect && (M==N))
+   {
+      real_t u = 0.0;
+      for (int i = 0; i < N; i++)
+      {
+         for (int j = i; j < N; j++)
+         {
+            real_t l = 0.0;
+            for (int k = 0; k < N; k++)
+            {
+               l += Z[i+k*N] * W[k] * Z[j+k*N];
+            }
+            l = fabs(l-data[i+N*j]);
+            if (l > u)
+            {
+               u = l;
+            }
+         }
+      }
+      if (u > 1e-9)
+      {
+         mfem::err << "dsyevr_Eigensystem(...):"
+                   << " max matrix deviation = " << u
+                   << "\n W[max] = " << W[N-1] << ", W[min] = "
+                   << W[0] << ", N = " << N << endl;
+      }
+      if (u > 1e-5)
+      {
+         mfem_error("dsyevr_Eigensystem(...): ERROR: ...");
+      }
+   }
+#endif
+   delete [] IWORK;
+   delete [] WORK;
+   delete [] ISUPPZ;
+   delete [] A;
 #else
-
+   MFEM_CONTRACT_VAR(a);
    MFEM_CONTRACT_VAR(ev);
    MFEM_CONTRACT_VAR(evect);
-   mfem_error("DenseMatrix::Eigensystem: Compiled without LAPACK");
-
 #endif
 }
 
+void dsyevr_Eigensystem(DenseMatrix &a, Vector &ev, DenseMatrix *evect)
+{
+  // dsyev_Eigensystem(a, ev, evect);
+   dsyevr_Eigensystem(a, ev, evect, 'A', 0.0, 1.0, 0, 1);
+}
+
+void dsyevr_Eigensystem(DenseMatrix &a, Vector &ev, DenseMatrix *evect,
+                        real_t VL, real_t VU)
+{
+   dsyevr_Eigensystem(a, ev, evect, 'V', VL, VU, 0, 1);
+}
+
+void dsyevr_Eigensystem(DenseMatrix &a, Vector &ev, DenseMatrix *evect,
+                        int IL, int IU)
+{
+   dsyevr_Eigensystem(a, ev, evect, 'I', 0.0, 1.0, IL, IU);
+}
+
 void dsygv_Eigensystem(DenseMatrix &a, DenseMatrix &b, Vector &ev,
                        DenseMatrix *evect)
 {
@@ -1185,14 +1220,14 @@ void dsygv_Eigensystem(DenseMatrix &a, DenseMatrix &b, Vector &ev,
       B[i] = b_data[i];
    }
 
-   MFEM_LAPACK_PREFIX(sygv_)(&ITYPE, &JOBZ, &UPLO, &N, A, &LDA, B, &LDB, W,
-                             &QWORK, &LWORK, &INFO);
+   MFEM_LAPACK_PREFIX(sygv_)(&ITYPE, &JOBZ, &UPLO, &N,
+      A, &LDA, B, &LDB, W, &QWORK, &LWORK, &INFO);
 
    LWORK = (int) QWORK;
    WORK = new real_t[LWORK];
 
-   MFEM_LAPACK_PREFIX(sygv_)(&ITYPE, &JOBZ, &UPLO, &N, A, &LDA, B, &LDB, W, WORK,
-                             &LWORK, &INFO);
+   MFEM_LAPACK_PREFIX(sygv_)(&ITYPE, &JOBZ, &UPLO, &N,
+      A, &LDA, B, &LDB, W, WORK, &LWORK, &INFO);
 
    if (INFO != 0)
    {
@@ -1211,13 +1246,202 @@ void dsygv_Eigensystem(DenseMatrix &a, DenseMatrix &b, Vector &ev,
 #endif
 }
 
+void dsygvx_Eigensystem(DenseMatrix &a, DenseMatrix &b,
+                        Vector &ev, DenseMatrix *evect,
+                        char RANGE, real_t VL, real_t VU, int IL, int IU)
+{
+#ifdef MFEM_USE_LAPACK
+   ev.SetSize(a.Width());
+   int       ITYPE    = 1;
+   char      JOBZ     = 'N';
+   char      UPLO     = 'U';
+   int       N        = a.Width();
+   real_t   *A        = new real_t[N*N];
+   int       LDA      = N;
+   real_t   *B        = new real_t[N*N];
+   int       LDB      = N;
+   real_t    ABSTOL   = 0.0;
+   int       M;
+   real_t   *W        = ev.GetData();
+   real_t   *Z        = NULL;
+   int       LDZ      = 1;
+   real_t    QWORK;
+   real_t   *WORK     = NULL;
+   int       LWORK    = -1; // query optimal (double) workspace size
+   int      *IWORK    = new int[5*N];
+   int      *IFAIL    = new int[N];
+   int       INFO;
+
+   if (evect) // Compute eigenvectors too
+   {
+      evect->SetSize(N);
+
+      JOBZ     = 'V';
+      Z        = evect->Data();
+      LDZ      = N;
+   }
+
+   int hw = a.Height() * a.Width();
+   real_t *data_a = a.Data();
+   real_t *data_b = b.Data();
+
+   for (int i = 0; i < hw; i++)
+   {
+      A[i] = data_a[i];
+      B[i] = data_b[i];
+   }
+
+   MFEM_LAPACK_PREFIX(sygvx_)( &ITYPE, &JOBZ, &RANGE, &UPLO,
+            &N, A, &LDA, B, &LDB, &VL, &VU, &IL, &IU,
+            &ABSTOL, &M, W, Z, &LDZ, &QWORK, &LWORK,
+            IWORK, IFAIL, &INFO );
+
+   LWORK  = (int) QWORK;
+   WORK  = new real_t[LWORK];
+
+   MFEM_LAPACK_PREFIX(sygvx_)( &ITYPE, &JOBZ, &RANGE, &UPLO,
+            &N, A, &LDA, B, &LDB, &VL, &VU, &IL, &IU,
+            &ABSTOL, &M, W, Z, &LDZ, WORK, &LWORK,
+            IWORK, IFAIL, &INFO );
+
+   if (INFO != 0)
+   {
+      mfem::err << "dsyevr_Eigensystem(...): DSYEVR error code: "
+                << INFO << endl;
+      mfem_error();
+   }
+
+   if (evect) // Compute eigenvectors too
+   {
+      evect->SetSize(N,M);
+   }
+   delete [] IFAIL;
+   delete [] IWORK;
+   delete [] WORK;
+   delete [] B;
+   if (evect == NULL) { delete [] A; }
+#else
+   MFEM_CONTRACT_VAR(a);
+   MFEM_CONTRACT_VAR(ev);
+   MFEM_CONTRACT_VAR(evect);
+#endif
+}
+
+void dsygvx_Eigensystem(DenseMatrix &a, DenseMatrix &b,
+                        Vector &ev, DenseMatrix *evect)
+{
+   //dsygv_Eigensystem(a, b, ev, evect);
+   dsygvx_Eigensystem(a, b, ev, evect, 'A', 0.0, 1.0, 0, 1);
+}
+
+void dsygvx_Eigensystem(DenseMatrix &a, DenseMatrix &b,
+                        Vector &ev, DenseMatrix *evect,
+                        real_t VL, real_t VU)
+{
+   dsygvx_Eigensystem(a, b, ev, evect, 'V', VL, VU, 0, 1);
+}
+
+void dsygvx_Eigensystem(DenseMatrix &a, DenseMatrix &b,
+                        Vector &ev, DenseMatrix *evect,
+                        int IL, int IU)
+{
+   dsygvx_Eigensystem(a, b, ev, evect, 'I', 0.0, 1.0, IL, IU);
+}
+
+real_t DenseMatrix::Eigenvalue(int i)
+{
+   MFEM_VERIFY(Height() == Width(), "a has to be a square matrix");
+
+   if (i < 0) { i = Width() - 1; }
+#ifdef MFEM_USE_LAPACK
+   Vector ev;
+   dsyevr_Eigensystem(*this, ev, NULL, i+1, i+1);
+   return ev[0];
+#else
+   if (i != Width() - 1)
+   {
+      mfem_error("DenseMatrix::NullSpace: Compiled without LAPACK");
+   }
+   // Use power method
+   int n = a.Width();
+   real_t alpha, eval_i = 0.0, eval_prev = 0.0;
+   int iter = 0;
+
+#ifdef MFEM_USE_SINGLE
+   const real_t rel_tol = 1e-7;
+#elif defined MFEM_USE_DOUBLE
+   const real_t rel_tol = 1e-14;
+#else
+   MFEM_ABORT("Floating point type undefined");
+#endif
+   Vector x_tmp(n), x(n);
+   x.Randomize(696383532);
+   do
+   {
+      a.Mult(x, x_tmp);
+      a.Mult(x_tmp, x);
+      eval_prev = eval_i;
+      eval_i = x.Norml2();
+      x *= 1.0/eval_i;
+      eval_i = sqrt(eval_i);
+      iter += 2;
+   }
+   while ((iter < 10000) && (fabs(eval_i - eval_prev)/fabs(eval_i) > rel_tol));
+   MFEM_VERIFY(fabs(eval_i - eval_prev)/fabs(eval_i) <= rel_tol,
+               "Inverse power method did not converge."
+               << "\n\t iter      = " << iter
+               << "\n\t eval_i    = " << eval_i
+               << "\n\t eval_prev = " << eval_prev
+               << "\n\t fabs(eval_i - eval_prev)/fabs(eval_i) = "
+               << fabs(eval_i - eval_prev)/fabs(eval_i));
+   return eval_i;
+#endif
+}
+
+real_t DenseMatrix::Eigenvalue(DenseMatrix &b, int i)
+{
+#ifdef MFEM_USE_LAPACK
+   if (i < 0) { i = Width() - 1; }
+   Vector ev;
+   dsygvx_Eigensystem(*this, b, ev, NULL, i+1, i+1);
+   return ev[0];
+#else
+   MFEM_CONTRACT_VAR(ns);
+   MFEM_CONTRACT_VAR(tol);
+   mfem_error("DenseMatrix::Eigenvalue: Compiled without LAPACK");
+   return 0.0;
+#endif
+}
+
+void DenseMatrix::NullSpace(DenseMatrix &ns, real_t tol)
+{
+#ifdef MFEM_USE_LAPACK
+   Vector ev;
+   dsyevr_Eigensystem(*this, ev, &ns, -tol, tol);
+#else
+   MFEM_CONTRACT_VAR(ns);
+   MFEM_CONTRACT_VAR(tol);
+   mfem_error("DenseMatrix::NullSpace: Compiled without LAPACK");
+#endif
+}
+
+void DenseMatrix::Eigensystem(Vector &ev, DenseMatrix *evect)
+{
+#ifdef MFEM_USE_LAPACK
+   dsyevr_Eigensystem(*this, ev, evect);
+#else
+   MFEM_CONTRACT_VAR(ev);
+   MFEM_CONTRACT_VAR(evect);
+   mfem_error("DenseMatrix::Eigensystem: Compiled without LAPACK");
+#endif
+}
+
 void DenseMatrix::Eigensystem(DenseMatrix &b, Vector &ev,
                               DenseMatrix *evect)
 {
 #ifdef MFEM_USE_LAPACK
-
-   dsygv_Eigensystem(*this, b, ev, evect);
-
+   dsygvx_Eigensystem(*this, b, ev, evect);
+   //dsygv_old_Eigensystem(*this, b, ev, evect);
 #else
    MFEM_CONTRACT_VAR(b);
    MFEM_CONTRACT_VAR(ev);
@@ -1477,6 +1701,21 @@ void DenseMatrix::Symmetrize()
    kernels::Symmetrize(Height(), Data());
 }
 
+bool DenseMatrix::IsSymmetric(real_t tol)
+{
+   for (int i = 0; i < Height(); i++)
+   {
+      real_t L = 0.0;
+      for (int j = i+1; j < Width(); j++)
+      {
+         if (fabs((*this)(i, j) - (*this)(j, i)) > tol) { return false ; }
+      }
+   }
+
+   return true;
+}
+
+
 void DenseMatrix::Lump()
 {
    for (int i = 0; i < Height(); i++)

linalg/densemat.hpp

@@ -272,35 +272,51 @@ public:
    /// Compute the square of the Frobenius norm of the matrix
    real_t FNorm2() const { real_t s, n2; FNorm(s, n2); return s*s*n2; }
 
-   /// Compute eigenvalues of A x = ev x where A = *this
+   /** Compute eigenvalues of A x = ev x where A = *this
+       A is assumed to be symmetric. */
    void Eigenvalues(Vector &ev)
    { Eigensystem(ev); }
 
-   /// Compute eigenvalues and eigenvectors of A x = ev x where A = *this
+   /** Compute ith eigenvalue of A x = ev x where A = *this
+       A is assumed to be symmetric. */
+   real_t Eigenvalue(int i = -1);
+
+   /** Compute eigenvalues and eigenvectors of A x = ev x where A = *this
+       A is assumed to be symmetric. */
    void Eigenvalues(Vector &ev, DenseMatrix &evect)
    { Eigensystem(ev, &evect); }
 
-   /// Compute eigenvalues and eigenvectors of A x = ev x where A = *this
+   /** Compute eigenvalues and eigenvectors of A x = ev x where A = *this
+       A is assumed to be symmetric. */
    void Eigensystem(Vector &ev, DenseMatrix &evect)
    { Eigensystem(ev, &evect); }
 
-   /** Compute generalized eigenvalues and eigenvectors of A x = ev B x,
-       where A = *this */
+   /** Compute generalized eigenvalues of A x = ev B x, where A = *this
+       A and B are assumed to be symmetric. */
    void Eigenvalues(DenseMatrix &b, Vector &ev)
    { Eigensystem(b, ev); }
 
-   /// Compute generalized eigenvalues of A x = ev B x, where A = *this
+   /** Compute ith eigenvalue of A x = ev B x where A = *this
+       A and B are assumed to be symmetric. */
+   real_t Eigenvalue(DenseMatrix &b, int i = -1);
+
+   /** Compute generalized eigenvalues and eigenvectors of A x = ev B x,
+       where A = *this. A and B are assumed to be symmetric.*/
    void Eigenvalues(DenseMatrix &b, Vector &ev, DenseMatrix &evect)
    { Eigensystem(b, ev, &evect); }
 
    /** Compute generalized eigenvalues and eigenvectors of A x = ev B x,
-       where A = *this */
+       where A = *this. A and B are assumed to be symmetric.*/
    void Eigensystem(DenseMatrix &b, Vector &ev, DenseMatrix &evect)
    { Eigensystem(b, ev, &evect); }
 
    void SingularValues(Vector &sv) const;
    int Rank(real_t tol) const;
 
+   /** Compute the Null Space of the matrix, such that A x = 0,
+       where A = *this* and x is a column of ns */
+   void NullSpace(DenseMatrix &ns, real_t tol);
+
    /// Return the i-th singular value (decreasing order) of NxN matrix, N=1,2,3.
    real_t CalcSingularvalue(const int i) const;
 
@@ -322,7 +338,6 @@ public:
    void SetCol(int c, const real_t* col);
    void SetCol(int c, const Vector &col);
 
-
    /// Set all entries of a row to the specified value.
    void SetRow(int row, real_t value);
    /// Set all entries of a column to the specified value.
@@ -346,6 +361,7 @@ public:
    void Transpose(const DenseMatrix &A);
    /// (*this) = 1/2 ((*this) + (*this)^t)
    void Symmetrize();
+   bool IsSymmetric(real_t tol = 1e-10);
 
    void Lump();
 

linalg/lapack.hpp

@@ -53,6 +53,13 @@ MFEM_LAPACK_PREFIX(sygv_) (int *ITYPE, char *JOBZ, char *UPLO, int * N,
                            real_t *A, int *LDA, real_t *B, int *LDB, real_t *W,
                            real_t *WORK, int *LWORK, int *INFO);
 extern "C" void
+MFEM_LAPACK_PREFIX(sygvx_)(int *ITYPE, char *JOBZ, char *RANGE, char *UPLO,
+                           int *N, double *A, int *LDA, double *B, int *LDB,
+                           double *VL, double *VU, int *IL, int *IU,
+                           double *ABSTOL, int *M, double *W, double *Z,
+                           int *LDZ, double *WORK, int *LWORK,int *IWORK,
+                           int *IFAIL, int *INFO);
+extern "C" void
 MFEM_LAPACK_PREFIX(gesvd_)(char *JOBU, char *JOBVT, int *M, int *N, real_t *A,
                            int *LDA, real_t *S, real_t *U, int *LDU, real_t *VT,
                            int *LDVT, real_t *WORK, int *LWORK, int *INFO);

miniapps/CMakeLists.txt

@@ -37,3 +37,4 @@ add_subdirectory(tribol)
 add_subdirectory(hooke)
 add_subdirectory(dpg)
 add_subdirectory(hdiv-linear-solver)
+add_subdirectory(stabilized)

miniapps/nurbs/CMakeLists.txt

@@ -45,6 +45,10 @@ add_mfem_miniapp(nurbs_solenoidal
   MAIN nurbs_solenoidal.cpp
   LIBRARIES mfem)
 
+add_mfem_miniapp(nurbs_biharm
+  MAIN nurbs_biharm.cpp
+  LIBRARIES mfem)
+
 if (MFEM_ENABLE_TESTING)
   add_test(NAME nurbs_ex1_1d_r1_o2_ser
     COMMAND $<TARGET_FILE:nurbs_ex1> -no-vis
@@ -199,6 +203,10 @@ if (MFEM_ENABLE_TESTING)
     COMMAND $<TARGET_FILE:nurbs_solenoidal> -no-vis
         -m ${PROJECT_SOURCE_DIR}/data/cube-nurbs.mesh -r 1 -o 2)
 
+  add_test(NAME nurbs_biharm_ser
+    COMMAND $<TARGET_FILE:nurbs_biharm> -no-vis
+        -m ${PROJECT_SOURCE_DIR}/data/square-nurbs.mesh -r 1 -o 2)
+
 endif()
 
 if (MFEM_USE_MPI)

miniapps/nurbs/nurbs_biharm.cpp

@@ -0,0 +1,408 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+//
+// Stabilized Convection-Diffusion
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include <list>
+
+
+using namespace std;
+using namespace mfem;
+
+real_t kappa_param = 1.0;
+
+real_t dif_fun(const Vector & x)
+{
+   return kappa_param;
+}
+
+real_t force_fun(const Vector & x)
+{
+   int d = x.Size();
+
+   real_t kappa = dif_fun(x);
+
+   real_t pi  = (real_t)(M_PI);
+
+   real_t sx = sin(pi*x[0]);
+   real_t sy = 1.0;
+   real_t sz = 1.0;
+
+   if (d >= 2)
+   {
+      sy = sin(pi*x[1]);
+   }
+   if (d >= 3)
+   {
+      sz = sin(pi*x[2]);
+   }
+
+   return d*d*kappa*pi*pi*pi*pi*sx*sy*sz;
+}
+
+real_t sol_fun(const Vector & x)
+{
+   real_t pi  = (real_t)(M_PI);
+
+   real_t sx = sin(pi*x[0]);
+   real_t sy = 1.0;
+   real_t sz = 1.0;
+
+   int d = x.Size();
+   if (d >= 2)
+   {
+      sy = sin(pi*x[1]);
+   }
+   if (d >= 3)
+   {
+      sz = sin(pi*x[2]);
+   }
+
+   return sx*sy*sz;
+}
+
+void grad_fun(const Vector & x, Vector & a)
+{
+   real_t pi  = (real_t)(M_PI);
+
+   real_t sx = sin(pi*x[0]);
+   real_t cx = cos(pi*x[0]);
+   real_t sy = 1.0;
+   real_t cy = 1.0;
+   real_t sz = 1.0;
+   real_t cz = 1.0;
+
+   int d = x.Size();
+   if (d >= 2)
+   {
+      sy = sin(pi*x[1]);
+      cy = cos(pi*x[1]);
+   }
+   if (d >= 3)
+   {
+      sz = sin(pi*x[2]);
+      cz = cos(pi*x[2]);
+   }
+
+   a[0] = pi*cx*sy;
+   a[1] = pi*sx*cy;
+}
+
+//----------------------------------------------------------
+real_t lap_fun(const Vector & x)
+{
+   real_t pi  = (real_t)(M_PI);
+
+   real_t sx = sin(pi*x[0]);
+   real_t cx = cos(pi*x[0]);
+   real_t sy = 1.0;
+   real_t cy = 1.0;
+   real_t sz = 1.0;
+   real_t cz = 1.0;
+
+   int d = x.Size();
+   if (d >= 2)
+   {
+      sy = sin(pi*x[1]);
+      cy = cos(pi*x[1]);
+   }
+   if (d >= 3)
+   {
+      sz = sin(pi*x[2]);
+      cz = cos(pi*x[2]);
+   }
+
+   return -d*pi*pi*sx*sy*sz;
+}
+
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../../data/square-nurbs.mesh";
+   const char *per_file  = "none";
+   const char *ref_file  = "";
+   int ref_levels = 0;
+   Array<int> master(0);
+   Array<int> slave(0);
+   bool static_cond = false;
+   bool visualization = false;
+   real_t penalty = -1;
+   Array<int> order(1);
+   order[0] = 2;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly, -1 for auto.");
+   args.AddOption(&per_file, "-p", "--per",
+                  "Periodic BCS file.");
+   args.AddOption(&ref_file, "-rf", "--ref-file",
+                  "File with refinement data");
+   args.AddOption(&master, "-pm", "--master",
+                  "Master boundaries for periodic BCs");
+   args.AddOption(&slave, "-ps", "--slave",
+                  "Slave boundaries for periodic BCs");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&kappa_param, "-k", "--kappa",
+                  "Sets the diffusion parameters, should be positive."
+                  " Negative values are replaced with function defined in source.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(mfem::out);
+      return 1;
+   }
+
+   args.PrintOptions(mfem::out);
+
+   if (order.Min()< 2)
+   {
+      mfem_error("Wrong order.");
+   }
+   // 2. Read the mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
+   //    the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 3. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement and knot insertion of knots defined
+   //    in a refinement file. We choose 'ref_levels' to be the largest number
+   //    that gives a final mesh with no more than 50,000 elements.
+   {
+      // Mesh refinement as defined in refinement file
+      if (mesh->NURBSext && (strlen(ref_file) != 0))
+      {
+         mesh->RefineNURBSFromFile(ref_file);
+      }
+
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh->UniformRefinement();
+      }
+      mesh->PrintInfo();
+   }
+
+   // 4. Define a finite element space on the mesh. Here we use continuous
+   //    Lagrange finite elements of the specified order. If order < 1, we
+   //    instead use an isoparametric/isogeometric space.
+   FiniteElementCollection *fec;
+   NURBSExtension *NURBSext = NULL;
+   int own_fec = 0;
+
+   if (mesh->NURBSext)
+   {
+      fec = new NURBSFECollection(order[0]);
+      own_fec = 1;
+
+      int nkv = mesh->NURBSext->GetNKV();
+      if (order.Size() == 1)
+      {
+         int tmp = order[0];
+         order.SetSize(nkv);
+         order = tmp;
+      }
+
+      if (order.Size() != nkv ) { mfem_error("Wrong number of orders set."); }
+      NURBSext = new NURBSExtension(mesh->NURBSext, order);
+
+      // Read periodic BCs from file
+      std::ifstream in;
+      in.open(per_file, std::ifstream::in);
+      if (in.is_open())
+      {
+         int psize;
+         in >> psize;
+         master.SetSize(psize);
+         slave.SetSize(psize);
+         master.Load(in, psize);
+         slave.Load(in, psize);
+         in.close();
+      }
+      master.Print();
+      slave.Print();
+      NURBSext->ConnectBoundaries(master,slave);
+   }
+   else if (order[0] == -1) // Isoparametric
+   {
+      if (mesh->GetNodes())
+      {
+         fec = mesh->GetNodes()->OwnFEC();
+         own_fec = 0;
+         mfem::out << "Using isoparametric FEs: " << fec->Name() << endl;
+      }
+      else
+      {
+         mfem::out <<"Mesh does not have FEs --> Assume order 1.\n";
+         fec = new H1_FECollection(1, dim);
+         own_fec = 1;
+      }
+   }
+   else
+   {
+      if (order.Size() > 1) { cout <<"Wrong number of orders set, needs one.\n"; }
+      fec = new H1_FECollection(abs(order[0]), dim);
+      own_fec = 1;
+   }
+
+   FiniteElementSpace *fespace = new FiniteElementSpace(mesh, NURBSext, fec);
+   mfem::out << "Number of finite element unknowns: "
+             << fespace->GetTrueVSize() << endl;
+
+   // 5. Determine the list of true (i.e. conforming) essential boundary dofs.
+   //    In this example, the boundary conditions are defined by marking all
+   //    the boundary attributes from the mesh as essential (Dirichlet) and
+   //    converting them to a list of true dofs.
+   Array<int> ess_tdof_list;
+   if (mesh->bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(mesh->bdr_attributes.Max());
+      ess_bdr = 1;
+
+      // Remove periodic BCs
+      for (int i = 0; i < master.Size(); i++)
+      {
+         ess_bdr[master[i]-1] = 0;
+         ess_bdr[slave[i]-1] = 0;
+      }
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // 6. Set up the linear form b(.) which corresponds to the right-hand side of
+   //    the FEM linear system, which in this case is (1,phi_i) where phi_i are
+   //    the basis functions in the finite element fespace.
+   ConstantCoefficient u_dir(0.0);
+
+   Coefficient *kappa_tmp;
+   if (kappa_param < 0.0)
+   {
+      kappa_tmp = new FunctionCoefficient(dif_fun);
+   }
+   else
+   {
+      kappa_tmp = new ConstantCoefficient(kappa_param);
+   }
+
+   Coefficient& kappa = *kappa_tmp;
+   FunctionCoefficient force(force_fun);
+
+   LinearForm *b = new LinearForm(fespace);
+   b->AddDomainIntegrator(new DomainLFIntegrator(force));
+   b->Assemble();
+
+   // 7. Define the solution vector x as a finite element grid function
+   //    corresponding to fespace. Initialize x with initial guess of zero,
+   //    which satisfies the boundary conditions.
+   GridFunction x(fespace);
+   x = 0.0;
+
+   // 8. Set up the bilinear form a(.,.) on the finite element space
+   //    corresponding to the Laplacian operator -Delta, by adding the Diffusion
+   //    domain integrator.
+   BilinearForm *a = new BilinearForm(fespace);
+   a->AddDomainIntegrator(new LaplaceLaplaceIntegrator(kappa));
+
+   // 9. Assemble the bilinear form and the corresponding linear system,
+   //    applying any necessary transformations such as: eliminating boundary
+   //    conditions, applying conforming constraints for non-conforming AMR,
+   //    static condensation, etc.
+   if (static_cond) { a->EnableStaticCondensation(); }
+   a->Assemble();
+
+   SparseMatrix A;
+   Vector B, X;
+   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
+
+   mfem::out << "Size of linear system: " << A.Height() << endl;
+
+#ifndef MFEM_USE_SUITESPARSE
+   // 10. Define a simple Jacobi preconditioner and use it to
+   //     solve the system A X = B with PCG.
+   GSSmoother M(A);
+   GMRES(A, M, B, X, 1, 2000, 2000, 1e-16, 0.0);
+#else
+   // 10. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
+   UMFPackSolver umf_solver;
+   umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
+   umf_solver.SetOperator(A);
+   umf_solver.Mult(B, X);
+#endif
+
+   // 11. Recover the solution as a finite element grid function.
+   a->RecoverFEMSolution(X, *b, x);
+
+   // 12. Save the refined mesh and the solution. This output can be viewed later
+   //     using GLVis: "glvis -m refined.mesh -g sol.gf".
+   {
+      ofstream mesh_ofs("refined.mesh");
+      mesh_ofs.precision(8);
+      mesh->Print(mesh_ofs);
+      ofstream sol_ofs("sol.gf");
+      sol_ofs.precision(8);
+      x.Save(sol_ofs);
+      sol_ofs.close();
+   }
+
+   // 13. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *mesh << x << flush;
+   }
+
+   // 14. Error computation
+   Vector norm(3);
+   int order_quad = 3*order.Max() + 4;
+   const IntegrationRule *irs[Geometry::NumGeom];
+   for (int i=0; i < Geometry::NumGeom; ++i)
+   {
+      irs[i] = &(IntRules.Get(i, order_quad));
+   }
+   FunctionCoefficient sol_cf(sol_fun);
+   VectorFunctionCoefficient grad_cf(mesh->Dimension(), grad_fun);
+   FunctionCoefficient lap_cf(lap_fun);
+
+   norm[0]= x.ComputeL2Error(sol_cf,irs);
+   norm[1]= x.ComputeGradError(&grad_cf, irs);
+   norm[2] = x.ComputeLaplaceError(&lap_cf, irs);
+
+   mfem::out << "|| x_h - x_ex ||            = " << norm[0]  << "\n";
+   mfem::out << "|| grad x_h - grad x_ex ||  = " << norm[1] << "\n";
+   mfem::out << "|| lap x_h - lap x_ex ||    = " << norm[2] << "\n";
+
+   // 15. Save data in the VisIt format
+   VisItDataCollection visit_dc("Biharm", mesh);
+   visit_dc.RegisterField("solution", &x);
+   visit_dc.Save();
+
+   // 16. Free the used memory.
+   delete fespace;
+   if (own_fec) { delete fec; }
+   delete mesh;
+
+   return 0;
+}
+

miniapps/nurbs/nurbs_ex1.cpp

@@ -53,96 +53,11 @@ public:
 inline bool operator==(const Data& d1,const Data& d2) { return (d1.x == d2.x); }
 inline bool operator <(const Data& d1,const Data& d2) { return (d1.x  < d2.x); }
 
-/** Class for integrating the bilinear form a(u,v) := (Q Laplace u, v) where Q
-    can be a scalar coefficient. */
-class Diffusion2Integrator: public BilinearFormIntegrator
-{
-private:
-#ifndef MFEM_THREAD_SAFE
-   Vector shape,laplace;
-#endif
-   Coefficient *Q;
-
-public:
-   /// Construct a diffusion integrator with coefficient Q = 1
-   Diffusion2Integrator() { Q = NULL; }
-
-   /// Construct a diffusion integrator with a scalar coefficient q
-   Diffusion2Integrator (Coefficient &q) : Q(&q) { }
-
-   /** Given a particular Finite Element
-       computes the element stiffness matrix elmat. */
-   virtual void AssembleElementMatrix(const FiniteElement &el,
-                                      ElementTransformation &Trans,
-                                      DenseMatrix &elmat)
-   {
-      int nd = el.GetDof();
-      int dim = el.GetDim();
-      real_t w;
-
-#ifdef MFEM_THREAD_SAFE
-      Vector shape(nd);
-      Vector laplace(nd);
-#else
-      shape.SetSize(nd);
-      laplace.SetSize(nd);
-#endif
-      elmat.SetSize(nd);
-
-      const IntegrationRule *ir = IntRule;
-      if (ir == NULL)
-      {
-         int order;
-         if (el.Space() == FunctionSpace::Pk)
-         {
-            order = 2*el.GetOrder() - 2;
-         }
-         else
-         {
-            order = 2*el.GetOrder() + dim - 1;
-         }
-
-         if (el.Space() == FunctionSpace::rQk)
-         {
-            ir = &RefinedIntRules.Get(el.GetGeomType(),order);
-         }
-         else
-         {
-            ir = &IntRules.Get(el.GetGeomType(),order);
-         }
-      }
-
-      elmat = 0.0;
-      for (int i = 0; i < ir->GetNPoints(); i++)
-      {
-         const IntegrationPoint &ip = ir->IntPoint(i);
-         Trans.SetIntPoint(&ip);
-         w = -ip.weight * Trans.Weight();
-
-         el.CalcShape(ip, shape);
-         el.CalcPhysLaplacian(Trans, laplace);
-
-         if (Q)
-         {
-            w *= Q->Eval(Trans, ip);
-         }
-
-         for (int jj = 0; jj < nd; jj++)
-         {
-            for (int ii = 0; ii < nd; ii++)
-            {
-               elmat(ii, jj) += w*shape(ii)*laplace(jj);
-            }
-         }
-      }
-   }
-
-};
 
 int main(int argc, char *argv[])
 {
    // 1. Parse command-line options.
-   const char *mesh_file = "../../data/star.mesh";
+   const char *mesh_file = "../../data/square-nurbs.mesh";
    const char *per_file  = "none";
    const char *ref_file  = "";
    int ref_levels = -1;
@@ -421,7 +336,7 @@ int main(int argc, char *argv[])
    }
    else
    {
-      a->AddDomainIntegrator(new Diffusion2Integrator(one));
+      a->AddDomainIntegrator(new LaplaceIntegrator(one, -1.0));
       a->AddBdrFaceIntegrator(new DGDiffusionIntegrator(mone, 0.0, 0.0), neu_bdr);
    }
 

miniapps/stabilized/CMakeLists.txt

@@ -0,0 +1,51 @@
+# Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+
+list(APPEND STAB_COMMON_SOURCES
+     stab_tau.cpp stab_condif.cpp stab_navsto.cpp)
+list(APPEND STAB_COMMON_HEADERS
+     stab_tau.hpp stab_condif.hpp stab_navsto.hpp manu.hpp skew.hpp)
+
+set(STAB_COMMON_FILES
+    EXTRA_SOURCES ${STAB_COMMON_SOURCES}
+    EXTRA_HEADERS ${STAB_COMMON_HEADERS})
+
+add_mfem_miniapp(condif
+    MAIN condif.cpp
+    ${STAB_COMMON_FILES}
+    LIBRARIES mfem)
+
+add_mfem_miniapp(navsto
+    MAIN navsto.cpp
+    ${STAB_COMMON_FILES}
+    LIBRARIES mfem)
+
+if (MFEM_ENABLE_TESTING)
+  add_test(NAME ex_condif
+    COMMAND $<TARGET_FILE:condif> -no-vis
+    -m ${PROJECT_SOURCE_DIR}/data/square-nurbs.mesh -r 1 -o 2)
+endif()
+
+if (MFEM_USE_MPI)
+  add_mfem_miniapp(navsto_p
+    MAIN navsto_p.cpp
+    ${STAB_COMMON_FILES}
+    LIBRARIES mfem)
+
+  if (MFEM_ENABLE_TESTING)
+    add_test(NAME navsto_np=4
+      COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
+      ${MPIEXEC_PREFLAGS} $<TARGET_FILE:navsto> -no-vis
+      ${MPIEXEC_POSTFLAGS}
+      -m ${PROJECT_SOURCE_DIR}/data/square-nurbs.mesh -r 1 -o 2)
+  endif()
+endif()

miniapps/stabilized/condif.cpp

@@ -0,0 +1,405 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+//
+// Stabilized Convection-Diffusion
+
+#include "stab_condif.hpp"
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include <list>
+
+using namespace std;
+using namespace mfem;
+
+real_t att_param = 1.0;
+real_t kappa_param = 1.0;
+real_t pi  = (real_t)(M_PI);
+
+using VectorFun = std::function<void(const Vector & x, Vector & a)>;
+using ScalarFun = std::function<real_t(const Vector & x)>;
+
+#include "skew.hpp"
+#include "manu.hpp"
+
+void evaluate1D(Vector &x, Vector &f, GridFunction *gf, int lod)
+{
+   // Get Mesh and Nodes gridfunction
+   Mesh *mesh = gf->FESpace()->GetMesh();
+   GridFunction *nodes = mesh->GetNodes();
+   if (!nodes)
+   {
+      nodes = new GridFunction(gf->FESpace());
+      mesh->GetNodes(*nodes);
+   }
+
+   // Evaluate
+   std::list<pair<real_t,real_t>> sol;
+   Vector vals,coords;
+   for (int i = 0; i <  mesh->GetNE(); i++)
+   {
+      int geom       = mesh->GetElementBaseGeometry(i);
+      RefinedGeometry *refined_geo = GlobGeometryRefiner.Refine(( Geometry::Type)geom, 1, lod);
+
+      gf->GetValues(i, refined_geo->RefPts, vals);
+      nodes->GetValues(i, refined_geo->RefPts, coords);
+
+      for (int j = 0; j < vals.Size(); j++)
+      {
+         sol.push_back(std::make_pair(coords[j],vals[j]));
+      }
+   }
+
+   // Sort and make unique
+   sol.sort();
+   sol.unique();
+
+   // Convert to Vectors
+   x.SetSize(sol.size());
+   f.SetSize(sol.size());
+   int i = 0;
+   for (std::list<pair<real_t,real_t>>::iterator d = sol.begin() ; d != sol.end(); ++d, i++)
+   {
+      x[i] = d->first;
+      f[i] = d->second;
+   }
+}
+
+StabType GetStabilisationType(int stype)
+{
+   switch (stype)
+   {
+      case GALERKIN:
+         mfem::out<<"Galerkin formulation"<<std::endl;
+         break;
+      case SUPG:
+         mfem::out<<"SUPG formulation"<<std::endl;
+         break;
+      case GLS:
+         mfem::out<<"GLS formulation"<<std::endl;
+         break;
+      case VMS:
+         mfem::out<<"VMS formulation"<<std::endl;
+         break;
+      default:
+         mfem::out<<"GAL"<<"\t"<<"SUPG"<<"\t"<<"GLS"<<"\t"<<"VMS"<<std::endl;
+         mfem::out<<GALERKIN<<"\t"<<SUPG<<"\t"<<GLS<<"\t"<<VMS<<std::endl;
+         mfem_error("Wrong formulation");
+   }
+   return (StabType) stype;
+}
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../../data/inline-quad.mesh";
+   const char *ref_file  = "";
+   int problem = 0;
+   int sstype = -2;
+   bool static_cond = false;
+   bool visualization = false;
+   int lod = 0;
+   real_t penalty = -1;
+   Array<int> order(1);
+   order[0] = 2;
+   int ref_levels = 0;
+
+   bool mono = true;
+
+   OptionsParser args(argc, argv);
+
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ref_file, "-rf", "--ref-file",
+                  "File with refinement data");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh.");
+   args.AddOption(&kappa_param , "-k", "--kappa",
+                  "Sets the diffusion parameters, should be positive.");
+   args.AddOption(&att_param , "-a", "--att",
+                  "Sets the velocity direction");
+   args.AddOption(&problem, "-p", "--problem",
+                  "Select the problem to solve:\n\t"
+                  "  0 = convection skew-to-the mesh\n\t"
+                  "  1 = manufactured solution\n");
+   args.AddOption(&sstype, "-s", "--stab", " Stabilization type:\n\t"
+                  " -2 = Galerkin\n\t"
+                  " -1 = GLS\n\t"
+                  "  0 = SUPG\n\t"
+                  "  1 = VMS\n");
+   args.AddOption(&mono, "-mo", "--mono", "-co",
+                  "--comp",
+                  "Use a monolithic integrator or a composed one.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&lod, "-lod", "--level-of-detail",
+                  "Refinement level for 1D solution output (0 means no output).");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+
+   }
+   args.PrintOptions(mfem::out);
+
+
+
+   // 2. Read the mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
+   //    the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 3. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement and knot insertion of knots defined
+   //    in a refinement file. We choose 'ref_levels' to be the largest number
+   //    that gives a final mesh with no more than 50,000 elements.
+   {
+      // Mesh refinement as defined in refinement file
+      if (mesh->NURBSext && (strlen(ref_file) != 0))
+      {
+         mesh->RefineNURBSFromFile(ref_file);
+      }
+
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh->UniformRefinement();
+      }
+      mesh->PrintInfo();
+   }
+
+   // 4. Define a finite element space on the mesh. Here we use continuous
+   //    Lagrange finite elements of the specified order. If order < 1, we
+   //    instead use an isoparametric/isogeometric space.
+   FiniteElementCollection *fec;
+   NURBSExtension *NURBSext = nullptr;
+   int own_fec = 1;
+
+   if (mesh->NURBSext)
+   {
+      fec = new NURBSFECollection(order[0]);
+
+      int nkv = mesh->NURBSext->GetNKV();
+      if (order.Size() == 1)
+      {
+         int tmp = order[0];
+         order.SetSize(nkv);
+         order = tmp;
+      }
+
+      if (order.Size() != nkv ) { mfem_error("Wrong number of orders set."); }
+      NURBSext = new NURBSExtension(mesh->NURBSext, order);
+   }
+   else if (order[0] == -1) // Isoparametric
+   {
+      if (mesh->GetNodes())
+      {
+         fec = mesh->GetNodes()->OwnFEC();
+         own_fec = 0;
+         mfem::out << "Using isoparametric FEs: " << fec->Name() << endl;
+      }
+      else
+      {
+         mfem::out <<"Mesh does not have FEs --> Assume order 1.\n";
+         fec = new H1_FECollection(1, dim);
+      }
+   }
+   else
+   {
+      if (order.Size() > 1) { cout <<"Wrong number of orders set, needs one.\n"; }
+      fec = new H1_FECollection(abs(order[0]), dim);
+   }
+
+   FiniteElementSpace *fespace = new FiniteElementSpace(mesh, NURBSext, fec);
+   mfem::out << "Number of finite element unknowns: "
+             << fespace->GetTrueVSize() << endl;
+
+   // 5. Determine the list of true (i.e. conforming) essential boundary dofs.
+   //    In this example, the boundary conditions are defined by marking all
+   //    the boundary attributes from the mesh as essential (Dirichlet) and
+   //    converting them to a list of true dofs.
+   Array<int> ess_tdof_list;
+   if (mesh->bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(mesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // 6. Set up the linear form b(.) which corresponds to the right-hand side of
+   //    the FEM linear system, which in this case is (1,phi_i) where phi_i are
+   //    the basis functions in the finite element fespace.
+   VectorFunctionCoefficient *adv, *grad;
+   FunctionCoefficient *kappa,*force, *sol, *lap;
+
+   if (problem == 0)
+   {
+      if (mesh->Dimension() != 2) mfem_error("Advection skew to the mesh needs a 2D mesh!");
+      adv = new VectorFunctionCoefficient(mesh->Dimension(), skew::adv);
+      kappa= new FunctionCoefficient(skew::kappa);
+
+      force = new FunctionCoefficient(skew::force);
+      sol = new FunctionCoefficient(skew::sol);
+      grad = new VectorFunctionCoefficient(mesh->Dimension(), skew::grad);
+      lap = new FunctionCoefficient(skew::laplace);
+
+   }
+   else if (problem == 1)
+   {
+      adv = new VectorFunctionCoefficient(mesh->Dimension(), manufactured::adv);
+      kappa= new FunctionCoefficient(manufactured::kappa);
+
+      force = new FunctionCoefficient(manufactured::force);
+      sol = new FunctionCoefficient(manufactured::sol);
+      grad = new VectorFunctionCoefficient(mesh->Dimension(), manufactured::grad);
+      lap = new FunctionCoefficient(manufactured::laplace);
+   }
+   else 
+   {
+      mfem_error("Incorrect problem!");
+   }
+
+   // 7. Define the solution vector x as a finite element grid function
+   //    corresponding to fespace. Initialize x with initial guess of zero,
+   //    which satisfies the boundary conditions.
+   GridFunction x(fespace);
+   x.ProjectCoefficient(*sol);
+
+   if (problem == 1)
+   {
+      Vector norm(3);
+      norm[0] = x.ComputeL2Error(*sol);
+      norm[1] = x.ComputeGradError(grad);
+      norm[2] = x.ComputeLaplaceError(lap);
+
+      mfem::out << "|| x_h - x_ex || = " << norm[0]  << "\n";
+      mfem::out << "|| grad x_h - grad x_ex ||    = " << norm[1] << "\n";
+      mfem::out << "|| lap x_h - lap x_ex ||      = " << norm[2] << "\n";
+   }
+
+   // 9. Assemble the bilinear form and the corresponding linear system,
+   //    applying any necessary transformations such as: eliminating boundary
+   //    conditions, applying conforming constraints for non-conforming AMR,
+   //    static condensation, etc.
+   StabType stype = GetStabilisationType(sstype);
+   FFH92Tau tau (adv, kappa, fespace);
+   StabConDifComposition stab_condif_comp(adv, kappa, force, &tau);
+
+   BilinearForm a(fespace);
+   LinearForm b(fespace);
+
+   if (mono)
+   {
+      a.AddDomainIntegrator(new StabConDifIntegrator(adv, kappa, force, &tau, stype));
+      b.AddDomainIntegrator(new StabConDifIntegrator(adv, kappa, force, &tau, stype));
+   }
+   else
+   {
+      stab_condif_comp.SetBilinearIntegrators(&a, stype);
+      stab_condif_comp.SetLinearIntegrators(&b, stype);
+   }
+
+   a.Assemble();
+   b.Assemble();
+
+   if (static_cond) { a.EnableStaticCondensation(); }
+   SparseMatrix A;
+   Vector B, X;
+
+   a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+
+   mfem::out << "Size of linear system: " << A.Height() << endl;
+
+#ifndef MFEM_USE_SUITESPARSE
+   // 10. Define a simple Jacobi preconditioner and use it to
+   //     solve the system A X = B with PCG.
+   GSSmoother M(A);
+   GMRES(A, M, B, X, 1, 2000, 2000, 1e-16, 0.0);
+#else
+   // 10. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
+   UMFPackSolver umf_solver;
+   umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
+   umf_solver.SetOperator(A);
+   umf_solver.Mult(B, X);
+#endif
+
+   // 11. Recover the solution as a finite element grid function.
+   a.RecoverFEMSolution(X, b, x);
+
+   // 12. Save the refined mesh and the solution. This output can be viewed later
+   //     using GLVis: "glvis -m refined.mesh -g sol.gf".
+   {
+      ofstream mesh_ofs("refined.mesh");
+      mesh_ofs.precision(8);
+      mesh->Print(mesh_ofs);
+      ofstream sol_ofs("sol.gf");
+      sol_ofs.precision(8);
+      x.Save(sol_ofs);
+      sol_ofs.close();
+   }
+
+   // 13. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *mesh << x << flush;
+   }
+
+   if (mesh->Dimension() == 1 && lod > 0)
+   {
+      Vector coord, val;
+      evaluate1D(coord, val, &x, lod);
+
+      ofstream sol_ofs("solution.dat");
+      for (int i = 0; i < x.Size();i++)
+      {
+         sol_ofs<<coord[i] <<"\t"<<val[i]<<endl;
+      }
+      sol_ofs.close();
+   }
+
+   // 14. Error computation
+   if (problem == 1)
+   {
+      Vector norm(3);
+      norm[0] = x.ComputeL2Error(*sol);
+      norm[1] = x.ComputeGradError(grad);
+      norm[2] = x.ComputeLaplaceError(lap);
+
+      mfem::out << "|| x_h - x_ex || = " << norm[0]  << "\n";
+      mfem::out << "|| grad x_h - grad x_ex ||    = " << norm[1] << "\n";
+      mfem::out << "|| lap x_h - lap x_ex ||      = " << norm[2] << "\n";
+   }
+
+   // 15. Save data in the VisIt format
+   VisItDataCollection visit_dc("condif", mesh);
+   visit_dc.RegisterField("solution", &x);
+   visit_dc.Save();
+
+   // 16. Free the used memory.
+   delete fespace;
+   if (own_fec) { delete fec; }
+   delete mesh;
+   delete adv, grad;
+   delete kappa, force, sol, lap;
+
+   return 0;
+}
+

miniapps/stabilized/manu.hpp

@@ -0,0 +1,114 @@
+namespace manufactured
+{
+
+//----------------------------------------------------------
+void adv(const Vector & x, Vector & a)
+{
+   a[1] = 1.0/(1.0 + att_param*att_param);
+   a[0] = sqrt(1.0 - a[1]*a[1]);
+}
+
+//----------------------------------------------------------
+real_t kappa(const Vector & x)
+{
+   return kappa_param;
+}
+
+//----------------------------------------------------------
+real_t force(const Vector & x)
+{
+   int d = x.Size();
+
+   Vector a(d);
+   adv(x, a);
+   real_t ax = a[0];
+   real_t ay = 0.0;
+   real_t az = 0.0;
+   real_t k = kappa(x);
+
+   real_t sx = sin(pi*x[0]);
+   real_t cx = cos(pi*x[0]);
+   real_t sy = 1.0;
+   real_t cy = 1.0;
+   real_t sz = 1.0;
+   real_t cz = 1.0;
+
+   if (d >= 2)
+   {
+      sy = sin(pi*x[1]);
+      cy = cos(pi*x[1]);
+      ay = a[1];
+   }
+   if (d >= 3)
+   {
+      sz = sin(pi*x[2]);
+      cz = cos(pi*x[2]);
+      az = a[2];
+   }
+
+   return ax*pi*cx*sy*sz
+        + ay*pi*sx*cy*sz
+        + az*pi*sx*sy*cz + d*k*pi*pi*sx*sy*sz;
+}
+
+//----------------------------------------------------------
+real_t sol(const Vector & x)
+{
+   real_t sx = sin(pi*x[0]);
+   real_t sy = 1.0;
+   real_t sz = 1.0;
+
+   int d = x.Size();
+   if (d >= 2)
+   {
+      sy = sin(pi*x[1]);
+   }
+   if (d >= 3)
+   {
+      sz = sin(pi*x[2]);
+   }
+
+   return sx*sy*sz;
+}
+
+
+//----------------------------------------------------------
+void grad(const Vector & x, Vector &grad)
+{
+   real_t sx = sin(pi*x[0]);
+   real_t sy = 1.0;
+   real_t sz = 1.0;
+
+   real_t gx = pi*cos(pi*x[0]);
+   real_t gy = 0.0;
+   real_t gz = 0.0;
+
+   grad[0] = gx;
+
+   int d = x.Size();
+   if (d >= 2)
+   {
+      sy = sin(pi*x[1]);
+      gy = pi*cos(pi*x[1]);
+
+      grad[0] = gx*sy;
+      grad[1] = sx*gy;
+   }
+   if (d >= 3)
+   {
+      sz = sin(pi*x[2]);
+      gz = pi*cos(pi*x[2]);
+
+      grad[0] = gx*sy*sz;
+      grad[1] = sx*gy*sz;
+      grad[2] = sx*sy*gz;
+   }
+}
+
+//----------------------------------------------------------
+real_t laplace(const Vector & x)
+{
+   return -x.Size()*pi*pi*sol(x);
+}
+
+}

miniapps/stabilized/navsto.cpp

@@ -0,0 +1,282 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+//
+// Stabilized  Navier-Stokes
+
+#include "stab_navsto.hpp"
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include <list>
+
+using namespace std;
+using namespace mfem;
+
+real_t kappa_param = 1.0;
+real_t pi  = (real_t)(M_PI);
+
+using VectorFun = std::function<void(const Vector & x, Vector & a)>;
+using ScalarFun = std::function<real_t(const Vector & x)>;
+
+void sol_fun(const Vector & x, Vector &sol)
+{
+   sol = 0.0;
+   if ((x[1] - 0.99 > 0.0) &&
+       (fabs(x[0] - 0.5) < 0.49) )
+   {
+      sol[0] = 1.0;
+   }
+}
+
+real_t kappa_fun(const Vector & x)
+{
+   return kappa_param;
+}
+
+void force_fun(const Vector & x, Vector &f)
+{
+   f = 0.0;
+  // f[0] = x[1]*(1.0-x[1])*x[0]*(1.0-x[0]);
+}
+
+StabType GetStabilisationType(int stype)
+{
+   switch (stype)
+   {
+      case GALERKIN:
+         mfem::out<<"Galerkin formulation"<<std::endl;
+         break;
+      case SUPG:
+         mfem::out<<"SUPG formulation"<<std::endl;
+         break;
+      case GLS:
+         mfem::out<<"GLS formulation"<<std::endl;
+         break;
+      case VMS:
+         mfem::out<<"VMS formulation"<<std::endl;
+         break;
+      default:
+         mfem::out<<"GAL"<<"\t"<<"SUPG"<<"\t"<<"GLS"<<"\t"<<"VMS"<<std::endl;
+         mfem::out<<GALERKIN<<"\t"<<SUPG<<"\t"<<GLS<<"\t"<<VMS<<std::endl;
+         mfem_error("Wrong formulation");
+   }
+   return (StabType) stype;
+}
+
+int main(int argc, char *argv[])
+{
+   // Parse command-line options.
+   const char *mesh_file = "../../data/inline-quad.mesh";
+   const char *ref_file  = "";
+   int problem = 0;
+   int sstype = -2;
+   bool static_cond = false;
+   bool visualization = false;
+
+   real_t penalty = -1;
+   int order = 1;
+   int ref_levels = 0;
+
+   bool mono = true;
+
+   OptionsParser args(argc, argv);
+
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ref_file, "-rf", "--ref-file",
+                  "File with refinement data");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order isoparametric space.");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh.");
+   args.AddOption(&kappa_param , "-k", "--kappa",
+                  "Sets the diffusion parameters, should be positive.");
+   args.AddOption(&problem, "-p", "--problem",
+                  "Select the problem to solve:\n\t"
+                  "  0 = convection skew-to-the mesh\n\t"
+                  "  1 = manufactured solution\n");
+   args.AddOption(&sstype, "-s", "--stab", " Stabilization type:\n\t"
+                  " -2 = Galerkin\n\t"
+                  " -1 = GLS\n\t"
+                  "  0 = SUPG\n\t"
+                  "  1 = VMS\n");
+   args.AddOption(&mono, "-mo", "--mono", "-co",
+                  "--comp",
+                  "Use a monolithic integrator or a composed one.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+
+   }
+   args.PrintOptions(mfem::out);
+
+   // Read the mesh from the given mesh file. We can handle triangular,
+   // quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
+   // the same code.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   // Refine the mesh to increase the resolution. In this example we do
+   // 'ref_levels' of uniform refinement and knot insertion of knots defined
+   // in a refinement file. We choose 'ref_levels' to be the largest number
+   // that gives a final mesh with no more than 50,000 elements.
+   {
+      // Mesh refinement as defined in refinement file
+      if (mesh.NURBSext && (strlen(ref_file) != 0))
+      {
+         mesh.RefineNURBSFromFile(ref_file);
+      }
+
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh.UniformRefinement();
+      }
+      mesh.PrintInfo();
+   }
+
+   // Define a finite element space on the mesh. Here we use continuous
+   // Lagrange finite elements of the specified order. If order < 1, we
+   // instead use an isoparametric/isogeometric space.
+   Array<FiniteElementCollection *> fecs(2);
+   fecs[0] = new H1_FECollection(order, dim);
+   fecs[1] = new H1_FECollection(order, dim);
+
+   Array<FiniteElementSpace *> spaces(2);
+   spaces[0] = new FiniteElementSpace(&mesh, fecs[0], dim);
+   spaces[1] = new FiniteElementSpace(&mesh, fecs[1]);
+
+   mfem::out << "Number of finite element unknowns:\n"
+             << "\tVelocity = "<<spaces[0]->GetTrueVSize() << endl
+             << "\tPressure = "<<spaces[1]->GetTrueVSize() << endl;
+   // Determine the list of true (i.e. conforming) essential boundary dofs.
+   // In this example, the boundary conditions are defined by marking all
+   // the boundary attributes from the mesh as essential (Dirichlet) and
+   // converting them to a list of true dofs.
+   Array<Array<int> *> ess_bdr(2);
+   Array<int> ess_tdof_list;
+
+   Array<int> ess_bdr_u(spaces[0]->GetMesh()->bdr_attributes.Max());
+   Array<int> ess_bdr_p(spaces[1]->GetMesh()->bdr_attributes.Max());
+
+   ess_bdr_p = 0;
+   ess_bdr_u = 1;
+
+   ess_bdr[0] = &ess_bdr_u;
+   ess_bdr[1] = &ess_bdr_p;
+
+   // Set up the linear form b(.) which corresponds to the right-hand side of
+   // the FEM linear system, which in this case is (1,phi_i) where phi_i are
+   // the basis functions in the finite element fespace.
+
+   // Define the solution vector xp as a finite element grid function
+   Array<int> bOffsets(3);
+   bOffsets[0] = 0;
+   bOffsets[1] = spaces[0]->GetTrueVSize();
+   bOffsets[2] = spaces[1]->GetTrueVSize();
+   bOffsets.PartialSum();
+
+   BlockVector xp(bOffsets);
+
+   GridFunction x_u(spaces[0]);
+   GridFunction x_p(spaces[1]);
+
+   x_u.MakeTRef(spaces[0], xp.GetBlock(0), 0);
+   x_p.MakeTRef(spaces[1], xp.GetBlock(1), 0);
+
+   VectorFunctionCoefficient sol(dim, sol_fun);
+
+   x_u.ProjectCoefficient(sol);
+   x_p = 0.0;
+
+   x_u.SetTrueVector();
+   x_p.SetTrueVector();
+
+   // Define the output
+   VisItDataCollection visit_dc("navsto", &mesh);
+   visit_dc.RegisterField("u", &x_u);
+   visit_dc.RegisterField("p", &x_p);
+   visit_dc.SetCycle(0);
+   visit_dc.Save();
+
+   // Define the problem parameters
+   FunctionCoefficient kappa(kappa_fun);
+   VectorFunctionCoefficient force(dim, force_fun);
+
+   // Define the stabilisation parameters
+   VectorGridFunctionCoefficient adv(&x_u);
+   ElasticInverseEstimateCoefficient invEst(spaces[0]);
+   FFH92Tau tau(&adv, &kappa, &invEst, 4.0);
+   FF91Delta delta(&adv, &kappa, &invEst);
+
+   tau.print = delta.print = true;
+
+   // Define the block nonlinear form
+   BlockNonlinearForm Hform(spaces);
+   Hform.AddDomainIntegrator(new StabInNavStoIntegrator(kappa, force, tau, delta));
+   Array<Vector *> rhs(2);
+   rhs = nullptr; // Set all entries in the array
+   Hform.SetEssentialBC(ess_bdr, rhs);
+
+   // Set up the preconditioner
+   JacobianPreconditioner jac_prec(bOffsets,
+                                   Array<Solver *>({new GSSmoother(0,5),
+                                                    new GSSmoother(0,5)}));
+
+   // Set up the Jacobian solver
+   GeneralResidualMonitor j_monitor("\t\t\t\tFGMRES", 25);
+   FGMRESSolver j_gmres;
+   j_gmres.iterative_mode = false;
+   j_gmres.SetRelTol(1e-2);
+   j_gmres.SetAbsTol(1e-12);
+   j_gmres.SetMaxIter(300);
+   j_gmres.SetPrintLevel(-1);
+   j_gmres.SetMonitor(j_monitor);
+   j_gmres.SetPreconditioner(jac_prec);
+
+   // Set up the newton solver
+   SystemResidualMonitor newton_monitor("Newton", 1, bOffsets, &visit_dc);
+   NewtonSolver newton_solver;
+   newton_solver.iterative_mode = true;
+   newton_solver.SetPrintLevel(-1);
+   newton_solver.SetMonitor(newton_monitor);
+   newton_solver.SetRelTol(1e-4);
+   newton_solver.SetAbsTol(1e-8);
+   newton_solver.SetMaxIter(25);
+   newton_solver.SetSolver(j_gmres);
+   newton_solver.SetOperator(Hform);
+
+   // Solve the Newton system
+   Vector zero;
+   newton_solver.Mult(zero, xp);
+
+   // Save data in the VisIt format
+   visit_dc.SetCycle(999999);
+   visit_dc.Save();
+
+   // Free the used memory.
+   for (int i = 0; i < fecs.Size(); ++i)
+   {
+      delete fecs[i];
+   }
+   for (int i = 0; i < spaces.Size(); ++i)
+   {
+      delete spaces[i];
+   }
+
+   return 0;
+}
+

miniapps/stabilized/navsto_p.cpp

@@ -0,0 +1,297 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+//
+// Stabilized  Navier-Stokes
+
+#include "stab_navsto.hpp"
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include <list>
+
+using namespace std;
+using namespace mfem;
+
+real_t kappa_param = 1.0;
+real_t pi  = (real_t)(M_PI);
+
+using VectorFun = std::function<void(const Vector & x, Vector & a)>;
+using ScalarFun = std::function<real_t(const Vector & x)>;
+
+void sol_fun(const Vector & x, Vector &sol)
+{
+   sol = 0.0;
+   if ((x[1] - 0.99 > 0.0) &&
+       (fabs(x[0] - 0.5) < 0.49) )
+   {
+      sol[0] = 1.0;
+   }
+}
+
+real_t kappa_fun(const Vector & x)
+{
+   return kappa_param;
+}
+
+void force_fun(const Vector & x, Vector &f)
+{
+   f = 0.0;
+  // f[0] = x[1]*(1.0-x[1])*x[0]*(1.0-x[0]);
+}
+
+StabType GetStabilisationType(int stype)
+{
+   switch (stype)
+   {
+      case GALERKIN:
+         mfem::out<<"Galerkin formulation"<<std::endl;
+         break;
+      case SUPG:
+         mfem::out<<"SUPG formulation"<<std::endl;
+         break;
+      case GLS:
+         mfem::out<<"GLS formulation"<<std::endl;
+         break;
+      case VMS:
+         mfem::out<<"VMS formulation"<<std::endl;
+         break;
+      default:
+         mfem::out<<"GAL"<<"\t"<<"SUPG"<<"\t"<<"GLS"<<"\t"<<"VMS"<<std::endl;
+         mfem::out<<GALERKIN<<"\t"<<SUPG<<"\t"<<GLS<<"\t"<<VMS<<std::endl;
+         mfem_error("Wrong formulation");
+   }
+   return (StabType) stype;
+}
+
+int main(int argc, char *argv[])
+{
+   // Initialize MPI and HYPRE.
+   Mpi::Init(argc, argv);
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();
+
+   // Parse command-line options.
+   const char *mesh_file = "../../data/inline-quad.mesh";
+   const char *ref_file  = "";
+   int problem = 0;
+   int sstype = -2;
+   bool static_cond = false;
+   bool visualization = false;
+
+   real_t penalty = -1;
+   int order = 1;
+   int ref_levels = 0;
+
+   bool mono = true;
+
+   OptionsParser args(argc, argv);
+
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ref_file, "-rf", "--ref-file",
+                  "File with refinement data");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order isoparametric space.");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh.");
+   args.AddOption(&kappa_param , "-k", "--kappa",
+                  "Sets the diffusion parameters, should be positive.");
+   args.AddOption(&problem, "-p", "--problem",
+                  "Select the problem to solve:\n\t"
+                  "  0 = convection skew-to-the mesh\n\t"
+                  "  1 = manufactured solution\n");
+   args.AddOption(&sstype, "-s", "--stab", " Stabilization type:\n\t"
+                  " -2 = Galerkin\n\t"
+                  " -1 = GLS\n\t"
+                  "  0 = SUPG\n\t"
+                  "  1 = VMS\n");
+   args.AddOption(&mono, "-mo", "--mono", "-co",
+                  "--comp",
+                  "Use a monolithic integrator or a composed one.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0) args.PrintUsage(cout);
+      return 1;
+   }
+   if (myid == 0) args.PrintOptions(cout);
+
+   // Read the mesh from the given mesh file. We can handle triangular,
+   // quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
+   // the same code.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   // Refine the mesh to increase the resolution. In this example we do
+   // 'ref_levels' of uniform refinement and knot insertion of knots defined
+   // in a refinement file. We choose 'ref_levels' to be the largest number
+   // that gives a final mesh with no more than 50,000 elements.
+   {
+      // Mesh refinement as defined in refinement file
+      if (mesh.NURBSext && (strlen(ref_file) != 0))
+      {
+         mesh.RefineNURBSFromFile(ref_file);
+      }
+
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh.UniformRefinement();
+      }
+      if (myid == 0) mesh.PrintInfo();
+   }
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+
+   // Define a finite element space on the mesh. Here we use continuous
+   // Lagrange finite elements of the specified order. If order < 1, we
+   // instead use an isoparametric/isogeometric space.
+   Array<FiniteElementCollection *> fecs(2);
+   fecs[0] = new H1_FECollection(order, dim);
+   fecs[1] = new H1_FECollection(order, dim);
+
+   Array<ParFiniteElementSpace *> spaces(2);
+   spaces[0] = new ParFiniteElementSpace(&pmesh, fecs[0], dim);//, Ordering::byVDIM);
+   spaces[1] = new ParFiniteElementSpace(&pmesh, fecs[1]);
+
+   Array<int> tdof(num_procs),udof(num_procs),pdof(num_procs);
+   tdof = 0;
+   tdof[myid] = spaces[0]->TrueVSize();
+   MPI_Reduce(tdof.GetData(), udof.GetData(), num_procs, MPI_INT, MPI_MAX, 0, MPI_COMM_WORLD);
+
+   tdof = 0;
+   tdof[myid] = spaces[1]->TrueVSize();
+   MPI_Reduce(tdof.GetData(), pdof.GetData(), num_procs, MPI_INT, MPI_MAX, 0, MPI_COMM_WORLD);
+
+   if (myid == 0)
+   {
+      mfem::out << "Number of finite element unknowns:\n";
+      mfem::out << "\tVelocity = "<<spaces[0]->GlobalTrueVSize() << endl;
+      mfem::out << "\tPressure = "<<spaces[1]->GlobalTrueVSize() << endl;
+      mfem::out << "Number of finite element unknowns per partition:\n";
+      mfem::out <<  "\tVelocity = ";udof.Print(mfem::out, num_procs);
+      mfem::out <<  "\tPressure = ";pdof.Print(mfem::out, num_procs);
+   }
+
+   // Mark all velocity boundary dofs as essential
+   Array<Array<int> *> ess_bdr(2);
+  // Array<int> ess_tdof_list;
+
+   Array<int> ess_bdr_u(spaces[0]->GetMesh()->bdr_attributes.Max());
+   Array<int> ess_bdr_p(spaces[1]->GetMesh()->bdr_attributes.Max());
+
+   ess_bdr_p = 0;
+   ess_bdr_u = 1;
+
+   ess_bdr[0] = &ess_bdr_u;
+   ess_bdr[1] = &ess_bdr_p;
+
+   // Define the solution vector xp as a finite element grid function
+   Array<int> bOffsets(3);
+   bOffsets[0] = 0;
+   bOffsets[1] = spaces[0]->TrueVSize();
+   bOffsets[2] = spaces[1]->TrueVSize();
+   bOffsets.PartialSum();
+
+   BlockVector xp(bOffsets);
+
+   ParGridFunction x_u(spaces[0]);
+   ParGridFunction x_p(spaces[1]);
+
+   VectorFunctionCoefficient sol(dim, sol_fun);
+   x_u.ProjectCoefficient(sol);
+   x_p = 0.0;
+
+   x_u.GetTrueDofs(xp.GetBlock(0));
+   x_p.GetTrueDofs(xp.GetBlock(1));
+
+   VisItDataCollection visit_dc("navsto", &pmesh);
+   visit_dc.RegisterField("u", &x_u);
+   visit_dc.RegisterField("p", &x_p);
+   visit_dc.SetCycle(0);
+   visit_dc.Save();
+
+   // Define the problem parameters
+   FunctionCoefficient kappa(kappa_fun);
+   VectorFunctionCoefficient force(dim, force_fun);
+
+   // Define the stabilisation parameters
+   VectorGridFunctionCoefficient adv(&x_u);
+   ElasticInverseEstimateCoefficient invEst(spaces[0]);
+   FFH92Tau tau(&adv, &kappa, &invEst, 4.0);
+   FF91Delta delta(&adv, &kappa, &invEst);
+
+   tau.print = delta.print = (myid == 0);
+
+   // Define the block nonlinear form
+   ParBlockNonlinearForm Hform(spaces);
+   Hform.AddDomainIntegrator(new StabInNavStoIntegrator(kappa, force, tau, delta));
+   Array<Vector *> rhs(2);
+   rhs = nullptr; // Set all entries in the array
+   Hform.SetEssentialBC(ess_bdr, rhs);
+
+   // Set up the preconditioner
+   JacobianPreconditioner jac_prec(bOffsets,
+                                   Array<Solver *>({new HypreSmoother(),
+                                                    new HypreSmoother()}));
+
+   // Set up the Jacobian solver
+   GeneralResidualMonitor j_monitor(MPI_COMM_WORLD,"\t\t\t\tFGMRES", 25);
+   FGMRESSolver j_gmres(MPI_COMM_WORLD);
+   j_gmres.iterative_mode = false;
+   j_gmres.SetRelTol(1e-2);
+   j_gmres.SetAbsTol(1e-12);
+   j_gmres.SetMaxIter(300);
+   j_gmres.SetPrintLevel(-1);
+   j_gmres.SetMonitor(j_monitor);
+   j_gmres.SetPreconditioner(jac_prec);
+
+   // Set up the newton solver
+   SystemResidualMonitor newton_monitor(MPI_COMM_WORLD,"Newton", 1, bOffsets, &visit_dc, &xp,
+                                        Array<ParGridFunction *>({&x_u, &x_p}));
+   NewtonSolver newton_solver(MPI_COMM_WORLD);
+   newton_solver.iterative_mode = true;
+   newton_solver.SetPrintLevel(-1);
+   newton_solver.SetMonitor(newton_monitor);
+   newton_solver.SetRelTol(1e-4);
+   newton_solver.SetAbsTol(1e-8);
+   newton_solver.SetMaxIter(25);
+   newton_solver.SetSolver(j_gmres);
+   newton_solver.SetOperator(Hform);
+
+   // Solve the Newton system
+   Vector zero;
+   newton_solver.Mult(zero, xp);
+
+   // Save data in the VisIt format
+   // Define the output
+   // Save data in the VisIt format
+   visit_dc.SetCycle(999999);
+   visit_dc.Save();
+
+   // Free the used memory.
+   for (int i = 0; i < fecs.Size(); ++i)
+   {
+      delete fecs[i];
+   }
+   for (int i = 0; i < spaces.Size(); ++i)
+   {
+      delete spaces[i];
+   }
+
+   return 0;
+}
+

miniapps/stabilized/skew.hpp

@@ -0,0 +1,46 @@
+namespace skew
+{
+
+//----------------------------------------------------------
+void adv(const Vector & x, Vector & a)
+{
+   a[1] = 1.0/(1.0 + att_param*att_param);
+   a[0] = sqrt(1.0 - a[1]*a[1]);
+}
+
+//----------------------------------------------------------
+real_t kappa(const Vector & x)
+{
+   return kappa_param;
+}
+
+//----------------------------------------------------------
+real_t force(const Vector & x)
+{
+   return 0.0;
+}
+
+//----------------------------------------------------------
+real_t sol(const Vector & x)
+{
+   if ((x[1] - x[0] - 0.2 < 0.0)
+      &(x[0] + x[1] -0.99 < 0.0))
+   {
+      return 1.0;
+   }
+   return 0.0;
+}
+
+//----------------------------------------------------------
+void grad(const Vector & x, Vector &grad)
+{
+   grad = 0.0;
+}
+
+//----------------------------------------------------------
+real_t laplace(const Vector & x)
+{
+   return 0.0;
+}
+
+}

miniapps/stabilized/stab_condif.cpp

@@ -0,0 +1,234 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "stab_condif.hpp"
+
+using namespace mfem;
+
+StabConDifIntegrator::StabConDifIntegrator(VectorCoefficient *a,
+                                           Coefficient *k,
+                                           Coefficient *f,
+                                           Tau *t, StabType s)
+ : adv(a), kappa(k), force(f), tau(t), stab(s), own_tau(false)
+{
+   if (tau == nullptr)
+   {
+      tau = new FFH92Tau(adv, kappa, 12.0);
+      own_tau = true;
+   }
+   else
+   {
+      tau->SetConvection(adv);
+      tau->SetDiffusion(kappa);
+   }
+}
+
+StabConDifIntegrator::~StabConDifIntegrator()
+{
+   if (own_tau) { delete tau; }
+}
+
+const IntegrationRule &StabConDifIntegrator::GetRule(
+   const FiniteElement &trial_fe,
+   const FiniteElement &test_fe,
+   ElementTransformation &Trans)
+{
+   int order = trial_fe.GetOrder() + test_fe.GetOrder();
+   return IntRules.Get(trial_fe.GetGeomType(), order);
+}
+
+void StabConDifIntegrator::AssembleElementMatrix(const FiniteElement &el,
+                                                 ElementTransformation &Trans,
+                                                 DenseMatrix &elmat )
+{
+   int nd = el.GetDof();
+   int dim = el.GetDim();
+   real_t w,k,t = 0;
+   Vector a(dim);
+
+   elmat.SetSize(nd);
+   shape.SetSize(nd);
+   dshape.SetSize(nd,dim);
+   adshape.SetSize(nd);
+   laplace.SetSize(nd);
+   trail.SetSize(nd);
+   test.SetSize(nd);
+
+   const IntegrationRule *ir = NonlinearFormIntegrator::IntRule ? NonlinearFormIntegrator::IntRule : &GetRule(el, el, Trans);
+
+   elmat = 0.0;
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      Trans.SetIntPoint (&ip);
+      w = Trans.Weight() * ip.weight;
+
+      // Calculate shapes
+      el.CalcPhysShape(Trans, shape);
+      el.CalcPhysDShape(Trans, dshape);
+
+      // Evaluate coefficients
+      k = kappa->Eval(Trans, ip);
+      adv->Eval(a, Trans, ip);
+
+      // Galerkin convection term
+      dshape.Mult(a, adshape);
+      AddMult_a_VWt(w, shape, adshape, elmat);
+
+      // Galerkin diffusion term
+      AddMult_a_AAt(w*k, dshape, elmat);
+
+      if (stab != GALERKIN)
+      {
+         // Calculate shapes
+         el.CalcPhysLaplacian(Trans, laplace);
+
+         // Evaluate coefficients
+         t = tau->Eval(Trans, ip);
+
+         // Stablization term
+         // - GLS:  stab = -1
+         // - SUPG: stab =  0
+         // - VMS:  stab = +1
+         add(adshape, stab*k, laplace, test);
+         add(adshape, -k, laplace, trail);
+         AddMult_a_VWt(w*t, test, trail, elmat);
+      }
+   }
+}
+
+void StabConDifIntegrator::AssembleRHSElementVect(const FiniteElement &el,
+                                      ElementTransformation &Trans,
+                                      Vector &elvect)
+{
+   int nd = el.GetDof();
+   int dim = el.GetDim();
+   real_t w,k,t,f;
+   Vector a(dim);
+
+   elvect.SetSize(nd);
+   shape.SetSize(nd);
+   dshape.SetSize(nd,dim);
+   adshape.SetSize(nd);
+   laplace.SetSize(nd);
+   test.SetSize(nd);
+
+   const IntegrationRule *ir = LinearFormIntegrator::IntRule ? LinearFormIntegrator::IntRule : &GetRule(el, el, Trans);
+
+   elvect = 0.0;
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      Trans.SetIntPoint (&ip);
+      w = Trans.Weight() * ip.weight;
+
+      // Calculate shapes
+      el.CalcPhysShape(Trans, shape);
+
+      // Evaluate coefficients
+      f = force->Eval(Trans, ip);
+
+      // Galerkin term
+      elvect.Add(w*f, shape);
+
+      if (stab != GALERKIN)
+      {
+         // Calculate shapes
+         el.CalcPhysDShape(Trans, dshape);
+         el.CalcPhysLaplacian(Trans, laplace);
+
+         // Evaluate coefficients
+         k = kappa->Eval(Trans, ip);
+         adv->Eval(a, Trans, ip);
+         t = tau->Eval(Trans, ip);
+
+         // Advective derivative
+         dshape.Mult(a, adshape);
+
+         // Stablization term
+         // - GLS:  stab = -1
+         // - SUPG: stab =  0
+         // - VMS:  stab = +1
+         add(adshape, stab*k, laplace, test);
+         elvect.Add(w*f*t, test);
+      }
+   }
+}
+
+
+StabConDifComposition::StabConDifComposition(VectorCoefficient *a,
+                                             Coefficient *k,
+                                             Coefficient *f,
+                                             Tau *t)
+  :  adv(a), kappa(k), force(f), tau(t), own_tau(false)
+{
+   if (tau == nullptr)
+   {
+      tau = new FFH92Tau(adv, kappa, 12.0);
+      own_tau = true;
+   }
+   else
+   {
+      tau->SetConvection(adv);
+      tau->SetDiffusion(kappa);
+   }
+
+   // SUPG coefficients
+   adv_tau = new ScalarVectorProductCoefficient(*tau, *adv);
+   adv_tau_force = new ScalarVectorProductCoefficient(*force, *adv_tau);
+   adv_tau_kappa = new ScalarVectorProductCoefficient(*kappa, *adv_tau);
+   adv_tau_adv = new OuterProductCoefficient(*adv_tau, *adv);
+
+   // GLS/VMS coefficients
+   kappa_tau = new ProductCoefficient(*kappa, *tau);
+   kappa_tau_kappa = new ProductCoefficient(*kappa_tau, *kappa);
+   kappa_tau_force = new ProductCoefficient(*kappa_tau, *force);
+}
+
+StabConDifComposition::~StabConDifComposition()
+{
+   if (own_tau) { delete tau; }
+   delete adv_tau, adv_tau_kappa, adv_tau_adv, adv_tau_force,
+          kappa_tau, kappa_tau_kappa,kappa_tau_force;
+}
+
+void StabConDifComposition::SetBilinearIntegrators(BilinearForm *a, StabType stype)
+{
+   a->AddDomainIntegrator(new ConservativeConvectionIntegrator(*adv));
+   a->AddDomainIntegrator(new DiffusionIntegrator(*kappa));
+   if (stype == GALERKIN) return;
+
+   // Add SUPG terms
+   a->AddDomainIntegrator(new DiffusionIntegrator(*adv_tau_adv));
+   a->AddDomainIntegrator(new GradLaplaceIntegrator(*adv_tau_kappa, 1.0));
+   if (stype == SUPG) return;
+
+   // Add VMS/GLS terms
+   real_t  s = (stype == GLS)? -1.0: 1.0;
+   a->AddDomainIntegrator(new LaplaceGradIntegrator(*adv_tau_kappa,-s));
+   a->AddDomainIntegrator(new LaplaceLaplaceIntegrator(*kappa_tau_kappa,-s));
+}
+
+void StabConDifComposition::SetLinearIntegrators(LinearForm *b, StabType stype)
+{
+   b->AddDomainIntegrator(new DomainLFIntegrator(*force));
+   if (stype == GALERKIN) return;
+
+   // Add SUPG terms
+   b->AddDomainIntegrator(new DomainLFGradIntegrator(*adv_tau_force));
+   if (stype == SUPG) return;
+
+   // Add VMS/GLS terms
+   real_t  s = (stype == GLS)? -1.0: 1.0;
+   b->AddDomainIntegrator(new DomainLFLaplaceIntegrator(*kappa_tau_force,-s));
+}
+
+

miniapps/stabilized/stab_condif.hpp

@@ -0,0 +1,130 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_STAB_CONDIF_HPP
+#define MFEM_STAB_CONDIF_HPP
+
+#include "mfem.hpp"
+#include "stab_tau.hpp"
+
+namespace mfem
+{
+
+/** This Class defines a monolithic integrator for stabilized multi-dimensional
+    convection-diffusion.
+
+     $(a \cdot \nabla u, v) + (\kappa \nabla u, \nabla v) 
+    + \sum (a \cdot \nabla u - \kappa \Delta u, \tau (a \cdot \nabla v + s \kappa \Delta v))_e$
+
+     $(f, \nabla v) 
+    + \sum (f, \tau (a \cdot \nabla v + s \kappa \Delta v))_e$
+*/
+class StabConDifIntegrator : public BilinearFormIntegrator,
+                             public LinearFormIntegrator
+{
+protected:
+   /// The advection field
+   VectorCoefficient *adv;
+   /// The diffusion parameter and force fields
+   Coefficient *kappa, *force;
+
+   /// The stabilization parameter
+   Tau *tau;
+   bool own_tau;
+
+   StabType stab;
+
+private:
+   Vector laplace, shape, adshape, trail, test;
+   DenseMatrix dshape;
+
+public:
+   StabConDifIntegrator(VectorCoefficient *a,
+                        Coefficient *k,
+                        Coefficient *f,
+                        Tau *t = nullptr, StabType s = GALERKIN);
+
+   ~StabConDifIntegrator();
+
+   virtual void AssembleElementMatrix(const FiniteElement &el,
+                                      ElementTransformation &Tr,
+                                      DenseMatrix &elmat);
+
+   static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
+                                         const FiniteElement &test_fe,
+                                         ElementTransformation &Trans);
+
+   virtual void AssembleRHSElementVect(const FiniteElement &el,
+                                       ElementTransformation &Tr,
+                                       Vector &elvect);
+
+   using LinearFormIntegrator::AssembleRHSElementVect;
+}; 
+
+/** This Class composes standard integrators to obtain a stabilized formulation for 
+    multi-dimensional convection-diffusion.
+
+     $(a \cdot \nabla u, v) + (\kappa \nabla u, \nabla v) 
+    + \sum (a \cdot \nabla u - \kappa \Delta u, \tau (a \cdot \nabla v + s \kappa \Delta v))_e$
+
+     $(f, \nabla v) 
+    + \sum (f, \tau (a \cdot \nabla v + s \kappa \Delta v))_e$
+*/
+class StabConDifComposition
+{
+
+private:
+   /// The advection field
+   VectorCoefficient *adv;
+   /// The diffusion parameter and force fields
+   Coefficient *kappa, *force;
+
+   /// The stabilization parameter
+   Tau *tau;
+   bool own_tau;
+
+   //// Helper coefficients for defining the weak forms
+   VectorCoefficient *adv_tau;
+   Coefficient *kappa_tau;
+
+   /// SUPG coefficients
+   VectorCoefficient *adv_tau_force;
+   VectorCoefficient *adv_tau_kappa;
+   MatrixCoefficient *adv_tau_adv;
+
+   /// GLS/VMS coefficients
+   Coefficient *kappa_tau_kappa;
+   Coefficient *kappa_tau_force;
+
+public:
+
+   /** Constructor
+       @a a: is the advection velocity field.
+       @a k: is the diffusion param field.
+       @a f: is the force field. */
+   StabConDifComposition(VectorCoefficient *a,
+                         Coefficient *k,
+                         Coefficient *f,
+                         Tau *t = nullptr);
+
+   /// Destructor
+   ~StabConDifComposition();
+
+   /// This method sets the integrators for the bilinearform
+   void SetBilinearIntegrators(BilinearForm *a, StabType s);
+
+   /// This method sets the integrators for the linearform
+   void SetLinearIntegrators(LinearForm *b, StabType s);
+};
+
+} // namespace mfem
+
+#endif

miniapps/stabilized/stab_navsto.cpp

@@ -0,0 +1,394 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "stab_navsto.hpp"
+
+using namespace mfem;
+
+StabInNavStoIntegrator::StabInNavStoIntegrator(Coefficient &mu_,
+                                               VectorCoefficient &force_,
+                                               Tau &t, Tau &d, StabType s)
+ : c_mu(&mu_), c_force(&force_), tau(&t), delta(&d), stab(s)
+{ }
+
+void StabInNavStoIntegrator::SetDim(int dim_)
+{
+  if (dim_ != dim)
+  {
+     dim = dim_;
+     u.SetSize(dim);
+     f.SetSize(dim);
+     res.SetSize(dim);
+     up.SetSize(dim);
+     grad_u.SetSize(dim);
+     hess_u.SetSize(dim, (dim*(dim+1))/2);
+     grad_p.SetSize(dim);
+     hmap.SetSize(dim,dim);
+
+     if (dim == 2)
+     {
+        hmap(0,0) = 0;
+        hmap(0,1) =  hmap(1,0) =  1;
+        hmap(1,1) = 2; 
+     }
+     else if (dim == 2)
+     {
+        hmap(0,0) = 0;
+        hmap(0,1) = hmap(1,0) = 1;
+        hmap(0,2) = hmap(2,0) = 2;
+        hmap(1,1) = 3;
+        hmap(1,2) = hmap(2,1) = 4;
+        hmap(2,2) = 5;
+     }
+     else
+     {
+        mfem_error("Only implemented for 2D and 3D");
+     }
+  }
+}
+
+real_t StabInNavStoIntegrator::GetElementEnergy(
+   const Array<const FiniteElement *>&el,
+   ElementTransformation &Tr,
+   const Array<const Vector *>&elfun)
+{
+   if (el.Size() != 2)
+   {
+      mfem_error("StabInNavStoIntegrator::GetElementEnergy"
+                 " has incorrect block finite element space size!");
+   }
+   SetDim(el[0]->GetDim());
+   int dof_u = el[0]->GetDof();
+
+   sh_u.SetSize(dof_u);
+   elf_u.UseExternalData(elfun[0]->GetData(), dof_u, dim);
+
+   int intorder = 2*el[0]->GetOrder();
+   const IntegrationRule &ir = IntRules.Get(el[0]->GetGeomType(), intorder);
+
+   real_t energy = 0.0;
+
+   for (int i = 0; i < ir.GetNPoints(); ++i)
+   {
+      const IntegrationPoint &ip = ir.IntPoint(i);
+      Tr.SetIntPoint(&ip);
+
+      real_t w = ip.weight * Tr.Weight();
+
+      el[0]->CalcPhysShape(Tr, sh_u);
+      elf_u.MultTranspose(sh_u, u);
+
+      energy += w*(u*u)/2;
+   }
+
+   return energy;
+}
+
+void StabInNavStoIntegrator::AssembleElementVector(
+   const Array<const FiniteElement *> &el,
+   ElementTransformation &Tr,
+   const Array<const Vector *> &elfun,
+   const Array<Vector *> &elvec)
+{
+   if (el.Size() != 2)
+   {
+      mfem_error("StabInNavStoIntegrator::AssembleElementVector"
+                 " has finite element space of incorrect block number");
+   }
+
+   int dof_u = el[0]->GetDof();
+   int dof_p = el[1]->GetDof();
+
+   SetDim(el[0]->GetDim());
+   int spaceDim = Tr.GetSpaceDim();
+   bool hess = (el[0]->GetDerivType() == (int) FiniteElement::HESS);
+   if (dim != spaceDim)
+   {
+      mfem_error("StabInNavStoIntegrator::AssembleElementVector"
+                 " is not defined on manifold meshes");
+   }
+   elvec[0]->SetSize(dof_u*dim);
+   elvec[1]->SetSize(dof_p);
+
+   *elvec[0] = 0.0;
+   *elvec[1] = 0.0;
+
+   elf_u.UseExternalData(elfun[0]->GetData(), dof_u, dim);
+   elv_u.UseExternalData(elvec[0]->GetData(), dof_u, dim);
+
+   sh_u.SetSize(dof_u);
+   shg_u.SetSize(dof_u, dim);
+   ushg_u.SetSize(dof_u);
+   shh_u.SetSize(dof_u, (dim*(dim+1))/2);
+   sh_p.SetSize(dof_p);
+   shg_p.SetSize(dof_p, dim);
+
+   int intorder = 2*el[0]->GetOrder();
+   const IntegrationRule &ir = IntRules.Get(el[0]->GetGeomType(), intorder);
+
+   for (int i = 0; i < ir.GetNPoints(); ++i)
+   {
+      const IntegrationPoint &ip = ir.IntPoint(i);
+      Tr.SetIntPoint(&ip);
+      real_t w = ip.weight * Tr.Weight();
+      real_t mu = c_mu->Eval(Tr, ip);
+      c_force->Eval(f, Tr, ip);
+
+      // Compute shape and interpolate
+      el[0]->CalcPhysShape(Tr, sh_u);
+      elf_u.MultTranspose(sh_u, u);
+
+      el[0]->CalcPhysDShape(Tr, shg_u);
+      shg_u.Mult(u, ushg_u);
+      MultAtB(elf_u, shg_u, grad_u);
+
+      if (hess)
+      {
+         el[0]->CalcPhysHessian(Tr,shh_u);
+         MultAtB(elf_u, shh_u, hess_u);
+      }
+      else
+      {
+         shh_u = 0.0;
+         hess_u = 0.0;
+      }
+
+      el[1]->CalcPhysShape(Tr, sh_p);
+      real_t p = sh_p*(*elfun[1]);
+
+      el[1]->CalcPhysDShape(Tr, shg_p);
+      shg_p.MultTranspose(*elfun[1], grad_p);
+
+      // Compute strong residual
+      grad_u.Mult(u,res);   // Add convection
+      res += grad_p;        // Add pressure
+      res -= f;             // Subtract force
+      for (int i = 0; i < dim; ++i)
+      {
+         for (int j = 0; j < dim; ++j)
+         {
+            res[j] -= mu*(hess_u(j,hmap(i,i)) + 
+                          hess_u(i,hmap(j,i))); // Add diffusion
+         }
+      }
+
+      // Compute stability params
+      real_t t = tau->Eval(Tr, ip);
+      real_t d = delta->Eval(Tr, ip);
+
+      // Compute momentum weak residual
+      flux.Diag(-p + d*grad_u.Trace(),dim);  // Add pressure & LSIC to flux
+      grad_u.Symmetrize();                   // Grad to strain
+      flux.Add(2*mu,grad_u);                 // Add stress to flux
+      AddMult_a_VVt(-1.0, u, flux);          // Add convection to flux
+      AddMult_a_VWt(t, res, u, flux);        // Add SUPG to flux --> check order u and res
+      AddMult_a_ABt(w, shg_u, flux, elv_u);  // Add flux term to rhs
+      AddMult_a_VWt(-w, sh_u, f,    elv_u);  // Add force term to rhs
+
+      // Compute momentum weak residual
+      elvec[1]->Add(w*grad_u.Trace(), sh_p); // Add Galerkin term
+      shg_p.Mult(res, sh_p);                 // PSPG help term
+      elvec[1]->Add(w*t, sh_p);              // Add PSPG term  - sign looks worng?
+   }
+}
+
+void StabInNavStoIntegrator::AssembleElementGrad(
+   const Array<const FiniteElement*> &el,
+   ElementTransformation &Tr,
+   const Array<const Vector *> &elfun,
+   const Array2D<DenseMatrix *> &elmats)
+{
+   int dof_u = el[0]->GetDof();
+   int dof_p = el[1]->GetDof();
+
+   SetDim(el[0]->GetDim());
+   bool hess = (el[0]->GetDerivType() == (int) FiniteElement::HESS);
+
+   elf_u.UseExternalData(elfun[0]->GetData(), dof_u, dim);
+
+   elmats(0,0)->SetSize(dof_u*dim, dof_u*dim);
+   elmats(0,1)->SetSize(dof_u*dim, dof_p);
+   elmats(1,0)->SetSize(dof_p, dof_u*dim);
+   elmats(1,1)->SetSize(dof_p, dof_p);
+
+   *elmats(0,0) = 0.0;
+   *elmats(0,1) = 0.0;
+   *elmats(1,0) = 0.0;
+   *elmats(1,1) = 0.0;
+
+   sh_u.SetSize(dof_u);
+   shg_u.SetSize(dof_u, dim);
+   ushg_u.SetSize(dof_u);
+   sh_p.SetSize(dof_p);
+   shg_p.SetSize(dof_p, dim);
+
+   int intorder = 2*el[0]->GetOrder();
+   const IntegrationRule &ir = IntRules.Get(el[0]->GetGeomType(), intorder);
+
+   for (int i = 0; i < ir.GetNPoints(); ++i)
+   {
+      const IntegrationPoint &ip = ir.IntPoint(i);
+      Tr.SetIntPoint(&ip);
+      real_t w = ip.weight * Tr.Weight();
+      real_t mu = c_mu->Eval(Tr, ip);
+      real_t t = tau->Eval(Tr, ip);
+      real_t d = delta->Eval(Tr, ip);
+
+      el[0]->CalcPhysShape(Tr, sh_u);
+      elf_u.MultTranspose(sh_u, u);
+
+      el[0]->CalcPhysDShape(Tr, shg_u);
+      MultAtB(elf_u, shg_u, grad_u);
+
+      shg_u.Mult(u, ushg_u);
+
+      el[1]->CalcPhysShape(Tr, sh_p);
+      real_t p = sh_p*(*elfun[1]);
+
+      el[1]->CalcPhysDShape(Tr, shg_p);
+      shg_p.MultTranspose(*elfun[1], grad_p);
+
+      // u,u block
+      for (int i_u = 0; i_u < dof_u; ++i_u)
+      {
+         for (int j_u = 0; j_u < dof_u; ++j_u)
+         {
+            // Diffusion
+            real_t mat = 0.0;
+            for (int dim_u = 0; dim_u < dim; ++dim_u)
+            {
+              mat += shg_u(i_u,dim_u)*shg_u(j_u,dim_u);
+            }
+            mat *= mu;
+
+            // Convection
+            mat -= ushg_u(i_u)*sh_u(j_u);      // Galerkin
+            mat += t*ushg_u(i_u)*ushg_u(j_u);  // SUPG
+
+            mat *= w;
+            for (int dim_u = 0; dim_u < dim; ++dim_u)
+            {
+               (*elmats(0,0))(i_u + dim_u*dof_u, j_u + dim_u*dof_u) += mat;
+            }
+
+            for (int i_dim = 0; i_dim < dim; ++i_dim)
+            {
+               for (int j_dim = 0; j_dim < dim; ++j_dim)
+               {
+                  (*elmats(0,0))(i_u + i_dim*dof_u, j_u + j_dim*dof_u) +=
+                     (mu + d)*shg_u(i_u,j_dim)*shg_u(j_u,i_dim)*w;
+               }
+            }
+         }
+      }
+
+      // u,p and p,u blocks
+      for (int i_p = 0; i_p < dof_p; ++i_p)
+      {
+         for (int j_u = 0; j_u < dof_u; ++j_u)
+         {
+            for (int dim_u = 0; dim_u < dim; ++dim_u)
+            {
+               (*elmats(0,1))(j_u + dof_u * dim_u, i_p) += (shg_p(i_p, dim_u)*t*ushg_u(j_u)
+                                                           -shg_u(j_u,dim_u)*sh_p(i_p))*w;
+               (*elmats(1,0))(i_p, j_u + dof_u * dim_u) +=  shg_u(j_u,dim_u)*sh_p(i_p)*w;
+            }
+         }
+      }
+
+      // p,p block
+      AddMult_a_AAt(w*t, shg_p, *elmats(1,1));
+   }
+}
+
+
+void GeneralResidualMonitor::MonitorResidual(int it, real_t norm,
+                                             const Vector &r, bool final)
+{
+   if (it == 0)
+   {
+      norm0 = norm;
+   }
+
+   if ((print_level > 0 &&  it%print_level == 0) || final)
+   {
+      mfem::out << prefix << " iteration " << std::setw(2) << it
+                << " : ||r|| = " << norm
+                << ",  ||r||/||r_0|| = " << 100*norm/norm0<<" % \n";
+   }
+}
+
+void SystemResidualMonitor::MonitorResidual(int it, real_t norm,
+                                            const Vector &r, bool final)
+{
+   if (dc && (it > 0))
+   {
+      if (rank > 1)
+      {
+         for (int i = 0; i < nvar; ++i)
+         {
+            pgf[i]->Distribute(xp->GetBlock(i));
+         }
+      }
+      dc->SetCycle(it);
+      dc->Save();
+   }
+
+   Vector vnorm(nvar);
+
+   for (int i = 0; i < nvar; ++i)
+   {
+       Vector r_i(r.GetData() + bOffsets[i], bOffsets[i+1] - bOffsets[i]);
+       if ( rank == 1 )
+       {
+          vnorm[i] = r_i.Norml2();
+       }
+       else
+       {
+          vnorm[i] = sqrt(InnerProduct(MPI_COMM_WORLD, r_i, r_i));
+       }
+       if (it == 0) norm0[i] = vnorm[i];
+   }
+
+   bool print = (print_level > 0 &&  it%print_level == 0) || final;
+   if (print)
+   {
+      mfem::out << prefix << " iteration " << std::setw(3) << it <<"\n"
+                << " ||r||  \t"<< "||r||/||r_0||  \n";
+      for (int i = 0; i < nvar; ++i)
+      {
+         mfem::out <<vnorm[i]<<"\t"<< 100*vnorm[i]/norm0[i]<<" % \n";
+      }
+   }
+}
+
+
+void JacobianPreconditioner::SetOperator(const Operator &op)
+{
+   BlockOperator *jacobian = (BlockOperator *) &op;
+
+   for (int i = 0; i < prec.Size(); ++i)
+   {
+      prec[i]->SetOperator(jacobian->GetBlock(i,i));
+      SetDiagonalBlock(i, prec[i]);
+   }
+
+   SetBlock(1,0, const_cast<Operator*>(&jacobian->GetBlock(1,0)));
+}
+
+JacobianPreconditioner::~JacobianPreconditioner()
+{
+   for (int i = 0; i < prec.Size(); ++i)
+   {
+      delete prec[i];
+   }
+}
+

miniapps/stabilized/stab_navsto.hpp

@@ -0,0 +1,242 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_STAB_NAVSTO_HPP
+#define MFEM_STAB_NAVSTO_HPP
+
+#include "mfem.hpp"
+#include "stab_tau.hpp"
+
+namespace mfem
+{
+
+/** Stabilized incompressible Navier-Stokes integrator
+    Start with Galerkin for stokes - done
+    Add convection - done
+    Modify diffusion - done
+    Add difussion to residual - done
+    CHECK NUMBERING HESSIAN -> NURBS = WRONG?? 2D = ok --> 3D???   --> DONE NEEDS CHECKING???
+    Inverse estimate check order -> done
+    Add force                    -> done
+    Parallel -->  done
+
+
+    Add supg   -  rhs done, jac conv + press --> ignore diffusion for now
+    Add pspg   -  rhs done, jac conv + press --> ignore diffusion for now
+    Add lsic    -  rhs done, jac conv + press --> ignore diffusion for now
+
+    Add correct inverse estimate -> done?? number does not coincide with H&C
+
+    Add VMS/GLS
+    Add selection option of different stab modes
+
+    Add Hessian check to inverse estimate
+
+
+    Check
+    - Hessian numbering in 3D
+    - Power method --> Laplack / null-space
+    - Elastic Inverse estimate
+
+    Leopoldo P. Franca, Sérgio L. Frey
+    Stabilized finite element methods:
+    II. The incompressible Navier-Stokes equations.
+    Computer Methods in Applied Mechanics and Engineering, 99(2-3), 209-233.
+
+    https://doi.org/10.1016/0045-7825(92)90041-H
+    https://www.sciencedirect.com/science/article/pii/004578259290041H
+
+*/
+class StabInNavStoIntegrator : public BlockNonlinearFormIntegrator
+{
+private:
+   Coefficient *c_mu;
+   VectorCoefficient *c_force;
+   Vector u, f, grad_p;
+   DenseMatrix flux;
+
+   DenseMatrix elf_u, elv_u;
+//   Vector elf_u, elv_u;//
+   DenseMatrix elf_p, elv_p;
+   Vector sh_u, ushg_u, sh_p;
+   DenseMatrix shg_u, shh_u, shg_p, grad_u, hess_u;
+   Array2D<int> hmap;
+
+   /// The stabilization parameters
+   StabType stab;
+   Tau *tau = nullptr;
+   Tau *delta = nullptr;
+   Vector res, up;
+
+   /// The advection field
+   VectorCoefficient *adv = nullptr; // tbd???
+
+   int dim = -1;
+   void SetDim(int dim);
+
+public:
+   StabInNavStoIntegrator(Coefficient &mu_,
+                          VectorCoefficient &force_,
+                          Tau &t, Tau &d,
+                          StabType s = GALERKIN);
+
+   virtual real_t GetElementEnergy(const Array<const FiniteElement *>&el,
+                                   ElementTransformation &Tr,
+                                   const Array<const Vector *> &elfun);
+
+   /// Perform the local action of the NonlinearFormIntegrator
+   virtual void AssembleElementVector(const Array<const FiniteElement *> &el,
+                                      ElementTransformation &Tr,
+                                      const Array<const Vector *> &elfun,
+                                      const Array<Vector *> &elvec);
+
+   /// Assemble the local gradient matrix
+   virtual void AssembleElementGrad(const Array<const FiniteElement*> &el,
+                                    ElementTransformation &Tr,
+                                    const Array<const Vector *> &elfun,
+                                    const Array2D<DenseMatrix *> &elmats);
+};
+
+class GeneralResidualMonitor : public IterativeSolverMonitor
+{
+public:
+   GeneralResidualMonitor(const std::string& prefix_, int print_lvl)
+      : prefix(prefix_)
+   {
+      print_level = print_lvl;
+      rank = 1;
+   }
+
+   GeneralResidualMonitor(MPI_Comm comm,
+                          const std::string& prefix_, int print_lvl)
+      : prefix(prefix_)
+   {
+#ifndef MFEM_USE_MPI
+      print_level = print_lvl;
+#else
+      MPI_Comm_rank(comm, &rank);
+      if (rank == 0)
+      {
+         print_level = print_lvl;
+      }
+      else
+      {
+         print_level = -1;
+      }
+#endif
+   }
+
+   virtual void MonitorResidual(int it, real_t norm, const Vector &r, bool final);
+
+private:
+   const std::string prefix;
+   int rank, print_level;
+   mutable real_t norm0;
+};
+
+class SystemResidualMonitor : public IterativeSolverMonitor
+{
+public:
+   SystemResidualMonitor(const std::string& prefix_,
+                         int print_lvl,
+                         Array<int> &offsets,
+                         DataCollection *dc_ = nullptr)
+      : prefix(prefix_), bOffsets(offsets), dc(dc_)
+   {
+      print_level = print_lvl;
+      nvar = bOffsets.Size()-1;
+      norm0.SetSize(nvar);
+      rank = 1;
+   }
+
+   SystemResidualMonitor(MPI_Comm comm,
+                         const std::string& prefix_,
+                         int print_lvl,
+                         Array<int> &offsets)
+      : prefix(prefix_), bOffsets(offsets), dc(nullptr), xp(nullptr)
+   {
+#ifndef MFEM_USE_MPI
+      print_level = print_lvl;
+      rank = 1;
+#else
+      MPI_Comm_rank(comm, &rank);
+      if (rank == 0)
+      {
+         print_level = print_lvl;
+      }
+      else
+      {
+         print_level = -1;
+      }
+#endif
+      nvar = bOffsets.Size()-1;
+      norm0.SetSize(nvar);
+   }
+   SystemResidualMonitor(MPI_Comm comm,
+                         const std::string& prefix_,
+                         int print_lvl,
+                         Array<int> &offsets,
+                         DataCollection *dc_,
+                         BlockVector *x,
+                         Array<ParGridFunction *> pgf_)
+      : prefix(prefix_), bOffsets(offsets), dc(dc_), xp(x), pgf(pgf_)
+   {
+#ifndef MFEM_USE_MPI
+      print_level = print_lvl;
+      rank = 1;
+#else
+      MPI_Comm_rank(comm, &rank);
+      if (rank == 0)
+      {
+         print_level = print_lvl;
+      }
+      else
+      {
+         print_level = -1;
+      }
+#endif
+      nvar = bOffsets.Size()-1;
+      norm0.SetSize(nvar);
+   }
+
+
+   virtual void MonitorResidual(int it, real_t norm, const Vector &r, bool final);
+
+private:
+   const std::string prefix;
+   int print_level, nvar, rank;
+   mutable Vector norm0;
+  // Offsets for extracting block vector segments
+   Array<int> &bOffsets;
+   DataCollection *dc;
+   BlockVector *xp;
+   Array<ParGridFunction *> pgf;
+};
+
+// Custom block preconditioner for the Jacobian
+class JacobianPreconditioner : public BlockLowerTriangularPreconditioner //BlockDiagonalPreconditioner
+{
+protected:
+   Array<Solver *> prec;
+public:
+   JacobianPreconditioner(Array<int> &offsets, Array<Solver *> p)
+     : BlockLowerTriangularPreconditioner (offsets), prec(p)
+   { MFEM_VERIFY(offsets.Size()-1 == p.Size(), ""); };
+
+   virtual void SetOperator(const Operator &op);
+
+   virtual ~JacobianPreconditioner();
+};
+
+
+} // namespace mfem
+
+#endif

miniapps/stabilized/stab_tau.cpp

@@ -0,0 +1,124 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "stab_tau.hpp"
+
+using namespace mfem;
+
+real_t FFH92Tau::GetElementSize(ElementTransformation &T)
+{
+   const DenseMatrix &dxdxi = T.Jacobian();
+   row.SetSize(dim);
+   h.SetSize(dim);
+
+   for (int i = 0; i < dim; i++)
+   {
+      dxdxi.GetRow(i, row);
+      h[i] = row.Norml2();
+   }
+   switch (dim)
+   {
+      case 1:
+         return h[0];
+      case 2:
+         return h[0]*h[1]*sqrt(2.0/(h[0]*h[0] + h[1]*h[1]));
+      case 3:
+         return h[0]*h[1]*h[2]*(3.0/(h[0]*h[0] + h[1]*h[1] + h[2]*h[2]));
+   }
+   mfem_error("Wrong dim!");
+   return -1.0;
+}
+
+real_t FFH92Tau::GetInverseEstimate(ElementTransformation &T,
+                                    const IntegrationPoint &ip, real_t scale)
+{
+   if (Ci>0.0)
+   {
+      return 1.0/Ci;
+   }
+   else
+   {
+      return 1.0/(invEst_cf->Eval(T,ip)*scale);
+   }
+}
+
+real_t FFH92Tau::Eval(ElementTransformation &T,
+                      const IntegrationPoint &ip)
+{
+   real_t k = kappa->Eval(T, ip);
+   adv->Eval(a, T, ip);
+   real_t hk = GetElementSize(T);
+   real_t ci = GetInverseEstimate(T, ip, hk*hk);
+   real_t mk = std::min(1.0/3.0, 2*ci);
+   real_t ap = a.Normlp(p);
+   // Prevent division by zero
+   ap = std::max(ap,std::numeric_limits<real_t>::min());
+   real_t pe = mk*ap*hk/(k_fac*k); // k_fac = 2 for CD and k_fac = 4 for NS
+   real_t xi = std::min(pe,1.0);
+   real_t tau = hk*xi/(2*ap);
+
+   if (print)
+   {
+      std::cout<<"\n==========================\n";
+      std::cout<<"  kappa = "<<k  <<std::endl;
+      std::cout<<"  adv   = "; a.Print(std::cout);
+      std::cout<<"  h     = "<<hk <<std::endl;
+      std::cout<<"  Ci    = "<<ci<<" "
+               <<( (Ci<0) ? "(Computed)" :"(Specified)")<<std::endl;
+      std::cout<<"  1/Ci  = "<<1.0/ci<<std::endl;
+      std::cout<<"  mk    = "<<mk <<std::endl;
+      std::cout<<"  |a|_p = "<<ap <<std::endl;
+      std::cout<<"  Pe    = "<<pe <<std::endl;
+      std::cout<<"  xi    = "<<xi <<std::endl;
+      std::cout<<"  tau   = "<<tau<<std::endl;
+      std::cout<<"  tau   = "<<mk*hk*hk/(8*k)<<std::endl;
+      std::cout<<"==========================\n\n";
+      print = false;
+   }
+
+   return tau;
+}
+
+real_t FF91Delta::Eval(ElementTransformation &T,
+                      const IntegrationPoint &ip)
+{
+   real_t k = kappa->Eval(T, ip);
+   adv->Eval(a, T, ip);
+   real_t hk = GetElementSize(T);
+   real_t ci = GetInverseEstimate(T, ip, hk*hk);
+   real_t mk = std::min(1.0/3.0, 2*ci);
+   real_t ap = a.Normlp(p); // Preventing division by zero not necessary
+   real_t pe = mk*ap*hk/(k_fac*k);
+   real_t xi = std::min(pe,1.0);
+   real_t delta = lambda*ap*hk*xi;
+
+   if (print)
+   {
+      std::cout<<"\n==========================\n";
+      std::cout<<"  kappa  = "<<k  <<std::endl;
+      std::cout<<"  adv    = "; a.Print(std::cout);
+      std::cout<<"  h      = "<<hk <<std::endl;
+      std::cout<<"  Ci     = "<<ci<<" "
+               <<( (Ci<0) ? "(Computed)" :"(Specified)")<<std::endl;
+      std::cout<<"  1/Ci   = "<<1.0/ci<<std::endl;
+      std::cout<<"  mk     = "<<mk <<std::endl;
+      std::cout<<"  |a|_p  = "<<ap <<std::endl;
+      std::cout<<"  Pe     = "<<pe <<std::endl;
+      std::cout<<"  xi     = "<<xi <<std::endl;
+      std::cout<<"  lambda = "<<lambda <<std::endl;
+      std::cout<<"  delta  = "<<delta<<std::endl;
+      std::cout<<"==========================\n\n";
+      print = false;
+   }
+
+   return delta;
+}
+

miniapps/stabilized/stab_tau.hpp

@@ -0,0 +1,298 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_STAB_TAU_HPP
+#define MFEM_STAB_TAU_HPP
+
+#include "mfem.hpp"
+
+namespace mfem
+{
+
+/// Enumerate to indicate the stabilisation type.
+enum StabType
+{
+   GALERKIN = -2,
+   SUPG = 0,
+   GLS = -1,
+   VMS = 1
+};
+
+/// This Class defines a generic stabilisation parameter.
+class Tau: public Coefficient
+{
+protected:
+   /// The advection field
+   VectorCoefficient *adv;
+   /// The diffusion parameter field
+   Coefficient *kappa;
+   /// Dimension of the problem
+   int dim;
+   /// Velocity vector
+   Vector a;
+
+public:
+   /** Construct a stabilized confection-diffusion integrator with:
+       - @a a the convection velocity
+       - @a d the diffusion coefficient*/
+   Tau (VectorCoefficient *a_, Coefficient *k) : adv(a_), kappa(k)
+   {
+      dim = adv->GetVDim();
+      a.SetSize(dim);
+   };
+   /// Simple Constructor
+   Tau () : adv(nullptr), kappa(nullptr) {};
+
+   /// Set the convection coefficient
+   virtual void SetConvection(VectorCoefficient *a_)
+   {
+      adv = a_;
+      dim = adv->GetVDim();
+      a.SetSize(dim);
+   };
+
+   /// Set the convection coefficient
+   virtual void SetDiffusion(Coefficient *k_) { kappa = k_; };
+
+   /// Flag for printing
+   bool print = true;
+};
+
+/** This Class defines the stabilisation parameter for the multi-dimensional
+    convection-diffusion problem.
+    When @a k_fac =2 the parameter is defined as given in:
+
+    Franca, L.P., Frey, S.L., & Hughes, T.J.R.
+    Stabilized finite element methods:
+    I. Application to the advective-diffusive model.
+    Computer Methods in Applied Mechanics and Engineering, 95(2), 253-276.
+
+    This also works for the convection-diffusion part of the navier-Stokes problem.
+    When @a k_fac = 4 the parameter is defined as given in:
+
+    Franca, L.P., Frey, S.L.,
+    Stabilized finite element methods:
+    II. The incompressible Navier-Stokes equations.
+    Computer Methods in Applied Mechanics and Engineering, 99(2-3), 209-233.
+*/
+class FFH92Tau: public Tau
+{
+protected:
+   /// User provided inverse estimate of the elements
+   real_t Ci = -1.0;
+
+   /// If @a Ci is negative the is inverse estimate computed
+   Coefficient *invEst_cf = nullptr;
+   bool own_ie = false;
+
+   // Routine to get the inverse estimate at each point
+   real_t GetInverseEstimate(ElementTransformation &T,
+                             const IntegrationPoint &ip, real_t scale = 1.0);
+
+   /// The norm used for the velocity vector
+   real_t p = 2.0;
+
+   /// The facor used for computing the element Peclet/Reynolds number
+   real_t k_fac = 2.0;
+
+   /// Temp variable
+   Vector row;
+
+   /// Element size in different directions
+   Vector h;
+
+    /** Returns element size according to:
+
+        Harari, I, & Hughes, T.J.R.
+        What are C and h?: Inequalities for the analysis and design of
+        finite element methods.
+        Computer methods in applied mechanics and engineering 97(2), 157-192.
+   */
+   real_t GetElementSize(ElementTransformation &T);
+
+public:
+   /** Construct a stabilized confection-diffusion integrator with:
+       - @a a the convection velocity
+       - @a d the diffusion coefficient
+       - @a ie_cf for computing the inverse estimates
+       - @a f factor for computing the element Pe/Re number (default = 2)
+       - @a p which norm to use for the velocity magnitude (default = 2) */
+    FFH92Tau (VectorCoefficient *a, Coefficient *k, Coefficient *ie_cf,
+              real_t f = 2.0, real_t norm_p = 2.0)
+      :  Tau(a,k), invEst_cf(ie_cf), Ci(-1.0), k_fac(f), p(norm_p) {};
+
+   /** Construct a stabilized confection-diffusion integrator with:
+       - @a a the convection velocity
+       - @a d the diffusion coefficient
+       - @a fes to provide to coefficient for computing the inverse estimates
+       - @a f factor for computing the element Pe/Re number (default = 2)
+       - @a p which norm to use for the velocity magnitude (default = 2) */
+    FFH92Tau (VectorCoefficient *a, Coefficient *k,
+              FiniteElementSpace *fes,
+              real_t f = 2.0, real_t norm_p = 2.0)
+      :  Tau(a,k), Ci(-1.0), k_fac(f), p(norm_p)
+   {
+      invEst_cf = new InverseEstimateCoefficient(fes);
+      own_ie = true;
+   };
+
+   /** Construct a stabilized confection-diffusion integrator with:
+       - @a a the convection velocity
+       - @a d the diffusion coefficient
+       - @a c_explicity provided inverse estimate (default = 1.0/12.0)
+       - @a f factor for computing the element Pe/Re number (default = 2)
+       - @a p which norm to use for the velocity magnitude (default = 2)*/
+    FFH92Tau (VectorCoefficient *a, Coefficient *k,
+              real_t c_ = 1.0/12.0, real_t f = 2.0, real_t norm_p = 2.0)
+      :  Tau(a,k), Ci(c_), k_fac(f), p(norm_p)
+   {
+      invEst_cf = NULL;
+   };
+
+   /** Construct a stabilized confection-diffusion integrator with:
+       - @a ie_cf for computing the inverse estimatestes
+       - @a f factor for Pe/Re definition (default = 2)
+       - @a p which norm to use for the velocity magnitude (default = 2)
+       Convection and diffusion need to be specified later using
+       SetConvection and SetDiffusion, respectivly*/
+   FFH92Tau (Coefficient *ie_cf,
+             real_t f = 2.0, real_t norm_p = 2.0)
+      : invEst_cf(ie_cf), Ci(-1.0), k_fac(f) , p(norm_p) {};
+
+   /** Construct a stabilized confection-diffusion integrator with:
+       - @a fes to provide to coefficient for computing the inverse estimates
+       - @a f factor for Pe/Re definition (default = 2)
+       - @a p which norm to use for the velocity magnitude (default = 2)
+       Convection and diffusion need to be specified later using
+       SetConvection and SetDiffusion, respectivly*/
+   FFH92Tau (FiniteElementSpace *fes,
+             real_t f = 2.0, real_t norm_p = 2.0)
+      : Ci(-1.0), k_fac(f) , p(norm_p)
+   {
+      invEst_cf = new InverseEstimateCoefficient(fes);
+      own_ie = true;
+   };
+
+   /** Construct a stabilized confection-diffusion integrator with:
+       - @a c_explicity provided inverse estimate (default = 1.0/12.0)
+       - @a f factor for computing the element Pe/Re number (default = 2)
+       - @a p which norm to use for the velocity magnitude (default = 2)
+       Convection and diffusion need to be specified later using
+       SetConvection and SetDiffusion, respectivly*/
+   FFH92Tau (real_t c_ = 1.0/12.0, real_t f = 2.0, real_t norm_p = 2.0)
+      : Ci(c_), k_fac(f), p(norm_p)
+   {
+      invEst_cf = NULL;
+   };
+
+   /// Evaluate the coefficient at @a ip.
+   virtual real_t Eval(ElementTransformation &T,
+                       const IntegrationPoint &ip) override;
+
+    // Destructor
+   ~FFH92Tau()
+   { if (own_ie) { delete invEst_cf; } }
+};
+
+/** This Class defines the stabilisation parameter for the multi-dimensional
+    convection-diffusion problem.
+
+    This also works for the convection-diffusion part of the navier-Stokes problem.
+
+    Franca, L.P., Frey, S.L.,
+    Stabilized finite element methods:
+    II. The incompressible Navier-Stokes equations.
+    Computer Methods in Applied Mechanics and Engineering, 99(2-3), 209-233.
+*/
+class FF91Delta: public FFH92Tau
+{
+protected:
+   /// Overall scalling parameter
+   real_t lambda = 1.0;
+
+public:
+   /** Construct a stabilized confection-diffusion integrator with:
+       - @a a the convection velocity
+       - @a d the diffusion coefficient
+       - @a ie_cf for computing the inverse estimates
+       - @a f factor for computing the element Pe/Re number (default = 2)
+       - @a p which norm to use for the velocity magnitude (default = 2) */
+    FF91Delta (VectorCoefficient *a, Coefficient *k, Coefficient *ie_cf,
+               real_t l = 1.0, real_t f = 2.0, real_t norm_p = 2.0)
+      :  FFH92Tau(a,k,ie_cf,f,norm_p), lambda(l){};
+
+   /** Construct a stabilized confection-diffusion integrator with:
+       - @a a the convection velocity
+       - @a d the diffusion coefficient
+       - @a fes for computing the inverse estimates
+       - @a l overall scalling factor for delta (default = 1)
+       - @a f factor for computing the element Pe/Re number (default = 4)
+       - @a p which norm to use for the velocity magnitude  (default = 2) */
+    FF91Delta (VectorCoefficient *a, Coefficient *k,
+               FiniteElementSpace *fes,
+               real_t l = 1.0,
+               real_t f = 4.0, real_t norm_p = 2.0)
+      :  FFH92Tau(a,k,fes,f,norm_p), lambda(l){};
+
+   /** Construct a stabilized confection-diffusion integrator with:
+       - @a a the convection velocity
+       - @a d the diffusion coefficient
+       - @a c_explicity provided inverse estimate (default = 1.0/12.0)
+       - @a l overall scalling factor for delta (default = 1)
+       - @a f factor for computing the element Pe/Re number (default = 4)
+       - @a p which norm to use for the velocity magnitude (default = 2)*/
+    FF91Delta (VectorCoefficient *a, Coefficient *k,
+               real_t c_ = 1.0/12.0, real_t l = 1.0,
+               real_t f = 4.0, real_t norm_p = 2.0)
+      : FFH92Tau(a,k,c_,f,norm_p), lambda(l){};
+      
+   /** Construct a stabilized confection-diffusion integrator with:
+       - @a ie_cf for computing the inverse estimatestes
+       - @a f factor for Pe/Re definition (default = 2)
+       - @a p which norm to use for the velocity magnitude (default = 2)
+       Convection and diffusion need to be specified later using
+       SetConvection and SetDiffusion, respectivly*/
+   FF91Delta (Coefficient *ie_cf,
+              real_t l = 1.0, real_t f = 2.0, real_t norm_p = 2.0)
+      : FFH92Tau(ie_cf,f,norm_p), lambda(l){};
+
+   /** Construct a stabilized confection-diffusion integrator with:
+       - @a fes for computing the inverse estimates
+       - @a f factor for Pe/Re definition (default = 2)
+       - @a p which norm to use for the velocity magnitude (default = 2)
+       Convection and diffusion need to be specified later using
+       SetConvection and SetDiffusion, respectivly*/
+   FF91Delta (FiniteElementSpace *fes,
+              real_t l = 1.0, real_t f = 4.0, real_t norm_p = 2.0)
+      : FFH92Tau(fes,f,norm_p), lambda(l){};
+
+   /** Construct a stabilized confection-diffusion integrator with:
+       - @a c_explicity provided inverse estimate (default = 1.0/12.0)
+       - @a l overall scalling factor for delta (default = 1)
+       - @a f factor for computing the element Pe/Re number (default = 4)
+       - @a p which norm to use for the velocity magnitude (default = 2)
+       Convection and diffusion need to be specified later using
+       SetConvection and SetDiffusion, respectivly*/
+   FF91Delta (real_t c_ = 1.0/12.0, real_t l = 1.0,
+              real_t f = 2.0, real_t norm_p = 2.0)
+      : FFH92Tau(c_,f,norm_p), lambda(l){};
+
+   /// Evaluate the coefficient at @a ip.
+   virtual real_t Eval(ElementTransformation &T,
+                       const IntegrationPoint &ip) override;
+
+    // Destructor
+   ~FF91Delta(){};
+};
+
+} // namespace mfem
+
+#endif

tests/unit/linalg/test_matrix_dense.cpp

@@ -686,7 +686,7 @@ TEST_CASE("Exponential", "[DenseMatrix]")
 
 #ifdef MFEM_USE_LAPACK
 
-enum class TestCase { GenEigSPD, GenEigGE, SVD};
+enum class TestCase { GenEigSPD, GenEigGE, SVD, NULLSPACE};
 std::string TestCaseName(TestCase testcase)
 {
    switch (testcase)
@@ -697,6 +697,8 @@ std::string TestCaseName(TestCase testcase)
          return "Generalized Eigenvalue problem for a general matrix";
       case TestCase::SVD:
          return "Singular Value Decomposition for a general matrix";
+      case TestCase::NULLSPACE:
+         return "NULL space for a general matrix";
    }
    return "";
 }
@@ -705,7 +707,7 @@ TEST_CASE("Eigensystem Problems",
           "[DenseMatrix]")
 {
    auto testcase = GENERATE(TestCase::GenEigSPD, TestCase::GenEigGE,
-                            TestCase::SVD);
+                            TestCase::SVD, TestCase::NULLSPACE);
 
    CAPTURE(TestCaseName(testcase));
 
@@ -840,6 +842,37 @@ TEST_CASE("Eigensystem Problems",
          REQUIRE(USVt.MaxMaxNorm() == MFEM_Approx(0.));
       }
       break;
+      case TestCase::NULLSPACE:
+      {
+
+         Vector ev, ev2;
+         DenseMatrix evect, evect2;
+         DenseMatrix A(M);
+         A.Symmetrize();
+
+         A.Eigenvalues( ev, evect);
+         mfem::out<<"ev = "; ev.Print(mfem::out,4);
+         int N = M.Width();
+
+         DenseMatrix ns;
+         int nss;
+         for (int i = 0; i < N; i++)
+         {
+            AddMult_a_VVt(-ev[i], Vector(evect.GetColumn(i),N), A);
+
+            A.Eigenvalues( ev2, evect2);
+            mfem::out<<"ev = "; ev2.Print(mfem::out,4);
+            A.NullSpace(ns, 1e-9);
+            mfem::out<<"null = "; ev2.Print(mfem::out,4);
+            REQUIRE(ns.Width() == i+1);
+            REQUIRE(A.Eigenvalue() - ev2(3) == MFEM_Approx(0.));
+            REQUIRE(A.Eigenvalue(0) - ev2(0) == MFEM_Approx(0.));
+            REQUIRE(A.Eigenvalue(1) - ev2(1) == MFEM_Approx(0.));
+            REQUIRE(A.Eigenvalue(2) - ev2(2) == MFEM_Approx(0.));
+            REQUIRE(A.Eigenvalue(3) - ev2(3) == MFEM_Approx(0.));
+         }
+      }
+      break;
    }
 }