small updates

.gitignore

@@ -411,6 +411,15 @@ miniapps/tribol/contact-patch-test
 miniapps/diag-smoothers/abs-l1-jacobi
 miniapps/diag-smoothers/mg-abs-l1-jacobi
 
+miniapps/hdg/hdg_advection
+miniapps/hdg/hdg_advectionp
+miniapps/hdg/hdg_poisson
+miniapps/hdg/hdg_poissonp
+miniapps/hdg/mesh*
+miniapps/hdg/sol*
+miniapps/hdg/*.mesh
+miniapps/hdg/*.gf
+
 # Unit test binary and outputs
 tests/unit/output_meshes
 tests/unit/unit_tests

fem/bilinearform.cpp

@@ -103,6 +103,11 @@ BilinearForm::BilinearForm (FiniteElementSpace * f, BilinearForm * bf, int ps)
    boundary_face_integs = bf->boundary_face_integs;
    boundary_face_integs_marker = bf->boundary_face_integs_marker;
 
+   /* HDG */
+   hdgintbfi = bf->hdgintbfi;
+   hdgbdrbfi = bf->hdgbdrbfi;
+
+
    AllocMat();
 }
 
@@ -273,6 +278,27 @@ void BilinearForm::AddBdrFaceIntegrator(BilinearFormIntegrator *bfi,
    boundary_face_integs_marker.Append(&bdr_marker);
 }
 
+/* HDG */
+void BilinearForm::AddHDGInteriorFaceIntegrator (BilinearFormIntegrator * bfi)
+{
+   hdgintbfi.Append (bfi);
+}
+
+/* HDG */
+void BilinearForm::AddHDGBoundaryFaceIntegrator (BilinearFormIntegrator * bfi)
+{
+   hdgbdrbfi.Append (bfi);
+   skeleton_boundary_face_integs_marker.Append(NULL);
+}
+
+/* HDG */
+void BilinearForm::AddHDGBoundaryFaceIntegrator (BilinearFormIntegrator * bfi,
+                                                 Array<int> &bdr_marker)
+{
+   hdgbdrbfi.Append (bfi);
+   skeleton_boundary_face_integs_marker.Append(&bdr_marker);
+}
+
 void BilinearForm::ComputeElementMatrix(int i, DenseMatrix &elmat) const
 {
    if (element_matrices)
@@ -749,6 +775,96 @@ void BilinearForm::Assemble(int skip_zeros)
       }
    }
 
+   /* HDG */
+   // Skeleton interior face integrals for HDG
+   if (hdgintbfi.Size())
+   {
+      FaceElementTransformations *ftr;
+      const FiniteElement *face_fe;
+      int nfaces = mesh->GetNumFaces();
+
+      // loop over all the edges
+      for (int i = 0; i < nfaces; i++)
+      {
+         ftr = mesh->GetInteriorFaceTransformations(i); // the transformation of the face
+         fes->GetFaceVDofs(i, vdofs);   // the defrees of freedom related to the face
+         face_fe = fes->GetFaceElement(
+                      i);   // point face_fe to the FiniteElement over the edge
+         if (ftr != NULL)
+         {
+            for (int k = 0; k < hdgintbfi.Size();
+                 k++) // Loop over the related interals, but there is only one hdgintbfi right now
+            {
+               hdgintbfi[k] -> AssembleFaceMatrix (*face_fe, *ftr,
+                                                   elemmat); // call AssembleFaceMatrix over the face
+               mat -> AddSubMatrix (vdofs, vdofs, elemmat,
+                                    skip_zeros);  // assemble the local matrix to the global one, skipping the zeros
+            }
+         }
+      }
+   }
+
+   /* HDG */
+   // Skeleton boundary face integrals for HDG
+   if (hdgbdrbfi.Size())
+   {
+      FaceElementTransformations *ftr;
+      const FiniteElement *face_fe;
+      // Which boundary attributes need to be processed?
+      Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
+                                 mesh->bdr_attributes.Max() : 0);
+      bdr_attr_marker = 0;
+      for (int k = 0; k < hdgbdrbfi.Size(); k++)
+      {
+         if (skeleton_boundary_face_integs_marker[k] == NULL)
+         {
+            bdr_attr_marker = 1;
+            break;
+         }
+         Array<int> &bdr_marker = *skeleton_boundary_face_integs_marker[k];
+         MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
+                     "invalid boundary marker for boundary face integrator #"
+                     << k << ", counting from zero");
+         for (int i = 0; i < bdr_attr_marker.Size(); i++)
+         {
+            bdr_attr_marker[i] |= bdr_marker[i];
+         }
+      }
+
+
+      int nbdrfaces = fes->GetNBE();
+
+      // loop over all the edges
+      for (int i = 0; i < nbdrfaces; i++)
+      {
+         const int bdr_attr = mesh->GetBdrAttribute(i);
+         if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
+
+         int face = mesh->GetBdrElementFaceIndex(i);
+         ftr = mesh->GetBdrFaceTransformations(i); // the transformation of the face
+
+         if (ftr != NULL)
+         {
+            fes->GetFaceVDofs(face, vdofs);   // the defrees of freedom related to the face
+            face_fe = fes->GetFaceElement(
+                         face);   // point face_fe to the FiniteElement over the edge
+            for (int k = 0; k < hdgbdrbfi.Size();
+                 k++) // Loop over the related interals, but there is only one hdgbdrbfi right now
+            {
+               if (skeleton_boundary_face_integs_marker[k] &&
+                   (*skeleton_boundary_face_integs_marker[k])[bdr_attr-1] == 0)
+               { continue; }
+
+               hdgbdrbfi[k] -> AssembleFaceMatrix (*face_fe, *ftr,
+                                                   elemmat); // call AssembleFaceMatrix over the face
+               mat -> AddSubMatrix (vdofs, vdofs, elemmat,
+                                    skip_zeros);  // assemble the local matrix to the global one, skipping the zeros
+            }
+         }
+      }
+   }
+
+
 #ifdef MFEM_USE_LEGACY_OPENMP
    if (free_element_matrices)
    {
@@ -1296,6 +1412,10 @@ BilinearForm::~BilinearForm()
       { delete interior_face_integs[k]; }
       for (k=0; k < boundary_face_integs.Size(); k++)
       { delete boundary_face_integs[k]; }
+      /* HDG */
+      // Delete skeleton integrals for HDG
+      for (k=0; k < hdgintbfi.Size(); k++) { delete hdgintbfi[k]; }
+      for (k=0; k < hdgbdrbfi.Size(); k++) { delete hdgbdrbfi[k]; }
    }
 }
 

fem/bilinearform.hpp

@@ -119,6 +119,15 @@ protected:
    Array<BilinearFormIntegrator*> boundary_face_integs;
    Array<Array<int>*> boundary_face_integs_marker; ///< Entries are not owned.
 
+   /* HDG */
+   /// Set of HDG skeleton face Integrators over interior face to be applied.
+   Array<BilinearFormIntegrator*> hdgintbfi;
+
+   /* HDG */
+   /// Set of HDG skeleton face Integrators over boundary face to be applied.
+   Array<BilinearFormIntegrator*> hdgbdrbfi;
+   Array<Array<int>*> skeleton_boundary_face_integs_marker;
+
    mutable DenseMatrix elemmat;
    mutable Array<int>  vdofs;
 
@@ -288,6 +297,14 @@ public:
    Array<Array<int>*> *GetBFBFI_Marker()
    { return &boundary_face_integs_marker; }
 
+   /* HDG */
+   // Array of the HDG type bilinear form integrators, right now there is only one
+   Array<BilinearFormIntegrator*> *GetHDGIntBFI() { return &hdgintbfi; }
+
+   /* HDG */
+   // Array of the HDG type bilinear form integrators, right now there is only one
+   Array<BilinearFormIntegrator*> *GetHDGBdrBFI() { return &hdgbdrbfi; }
+
    /// Returns a reference to: $ M_{ij} $
    const real_t &operator()(int i, int j) { return (*mat)(i,j); }
 
@@ -437,6 +454,17 @@ public:
    void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi,
                              Array<int> &bdr_marker);
 
+   /* HDG */
+   /// Adds HDG Interior Integrator.
+   void AddHDGInteriorFaceIntegrator(BilinearFormIntegrator *bfi);
+
+   /* HDG */
+   /// Adds new HDG Boundary Integrator.
+   void AddHDGBoundaryFaceIntegrator(BilinearFormIntegrator *bfi);
+
+   void AddHDGBoundaryFaceIntegrator(BilinearFormIntegrator *bfi,
+                                     Array<int> &bdr_marker);
+
    /// Sets all sparse values of $ M $ and $ M_e $ to 'a'.
    void operator=(const real_t a)
    {

fem/bilininteg.cpp

@@ -220,6 +220,76 @@ void BilinearFormIntegrator::AssembleFaceMatrix(
               " Integrator class.");
 }
 
+/* HDG optimized integrators starts */
+void BilinearFormIntegrator::AssembleElementMatrix2FES(const FiniteElement
+                                                       &fe_q,
+                                                       const FiniteElement &fe_u,
+                                                       ElementTransformation &Trans,
+                                                       DenseMatrix &elmat1)
+{
+   MFEM_ABORT("AssembleElementMatrix2FES is not implemented for this"
+              " Integrator class.");
+}
+
+void BilinearFormIntegrator::AssembleFaceMatrixOneElement1and1FES(
+   const FiniteElement &fe_u,
+   const FiniteElement &face_fe,
+   FaceElementTransformations &Trans,
+   const int elem1or2,
+   const bool onlyB,
+   DenseMatrix &elmat1,
+   DenseMatrix &elmat2,
+   DenseMatrix &elmat3,
+   DenseMatrix &elmat4)
+{
+   MFEM_ABORT("AssembleFaceMatrixOneElement1and1FES is not implemented for this"
+              " Integrator class.");
+}
+
+void BilinearFormIntegrator::AssembleFaceMatrixOneElement2and1FES(
+   const FiniteElement &fe_q,
+   const FiniteElement &fe_u,
+   const FiniteElement &face_fe,
+   FaceElementTransformations &Trans,
+   const int elem1or2,
+   const bool onlyB,
+   DenseMatrix &elmat1,
+   DenseMatrix &elmat2,
+   DenseMatrix &elmat3,
+   DenseMatrix &elmat4)
+{
+   MFEM_ABORT("AssembleFaceMatrixOneElement2and1FES is not implemented for this"
+              " Integrator class.");
+}
+
+void BilinearFormIntegrator::AssembleFaceMatrixOneElement2and2FES(
+   const FiniteElement &fe_q,
+   const FiniteElement &fe_u,
+   const FiniteElement &face_fe,
+   const FiniteElement &face_fe2,
+   FaceElementTransformations &Trans,
+   const int elem1or2,
+   const bool onlyB,
+   DenseMatrix &elmat1,
+   DenseMatrix &elmat2,
+   DenseMatrix &elmat3,
+   DenseMatrix &elmat4)
+{
+   MFEM_ABORT("AssembleFaceMatrixOneElement2and2FES is not implemented for this"
+              " Integrator class.");
+}
+/* HDG optimized integrators ends */
+
+/* HDG */
+void BilinearFormIntegrator::AssembleFaceMatrix(const FiniteElement &face_fe,
+                                                FaceElementTransformations &Trans,
+                                                DenseMatrix &elmat)
+{
+   MFEM_ABORT("AssembleFaceMatrix(const FiniteElement &face_fe,"
+              " FaceElementTransformations &Trans,"
+              " DenseMatrix &elmat) is not implemented for this");
+}
+
 void BilinearFormIntegrator::AssembleTraceFaceMatrix (int elem,
                                                       const FiniteElement &trial_face_fe,
                                                       const FiniteElement &test_fe1,

fem/bilininteg.hpp

@@ -194,6 +194,54 @@ public:
                                    FaceElementTransformations &Trans,
                                    DenseMatrix &elmat);
 
+   /* HDG */
+   /* For the case when there are 2 finite element spaces - such as LHDG */
+   virtual void AssembleElementMatrix2FES(const FiniteElement &fe_q,
+                                          const FiniteElement &fe_u,
+                                          ElementTransformation &Trans,
+                                          DenseMatrix &elmat1);
+
+   /* HDG */
+   /* For the optimized HDG calculations  */
+   /* 1 element based 1 face based FES */
+   virtual void AssembleFaceMatrixOneElement1and1FES(const FiniteElement &fe_u,
+                                                     const FiniteElement &face_fe,
+                                                     FaceElementTransformations &Trans,
+                                                     const int elem1or2,
+                                                     const bool onlyB,
+                                                     DenseMatrix &elmat1,
+                                                     DenseMatrix &elmat2,
+                                                     DenseMatrix &elmat3,
+                                                     DenseMatrix &elmat4);
+   /* 2 element based 1 face based FES */
+   virtual void AssembleFaceMatrixOneElement2and1FES(const FiniteElement &fe_q,
+                                                     const FiniteElement &fe_u,
+                                                     const FiniteElement &face_fe,
+                                                     FaceElementTransformations &Trans,
+                                                     const int elem1or2,
+                                                     const bool onlyB,
+                                                     DenseMatrix &elmat1,
+                                                     DenseMatrix &elmat2,
+                                                     DenseMatrix &elmat3,
+                                                     DenseMatrix &elmat4);
+   /* 2 element based 2 face based FES */
+   virtual void AssembleFaceMatrixOneElement2and2FES(const FiniteElement &fe_q,
+                                                     const FiniteElement &fe_u,
+                                                     const FiniteElement &face_fe,
+                                                     const FiniteElement &face_fe2,
+                                                     FaceElementTransformations &Trans,
+                                                     const int elem1or2,
+                                                     const bool onlyB,
+                                                     DenseMatrix &elmat1,
+                                                     DenseMatrix &elmat2,
+                                                     DenseMatrix &elmat3,
+                                                     DenseMatrix &elmat4);
+   // Assemble a local matrix over an edge, HDG skeleton integral
+   virtual void AssembleFaceMatrix(const FiniteElement &face_fe,
+                                   FaceElementTransformations &Trans,
+                                   DenseMatrix &elmat);
+   /* End of HDG integrators */
+
    /** Abstract method used for assembling TraceFaceIntegrators for
        DPG weak formulations. */
    virtual void AssembleTraceFaceMatrix(int elem,

fem/coefficient.cpp

@@ -2060,4 +2060,42 @@ CoefficientVector::~CoefficientVector()
    delete qf;
 }
 
+/* HDG */
+double ComputeMean(Coefficient &coeff, Mesh &mesh,
+                   const IntegrationRule *irs[])
+{
+   double norm = 0.0;
+   ElementTransformation *tr;
+
+   for (int i = 0; i < mesh.GetNE(); i++)
+   {
+      tr = mesh.GetElementTransformation(i);
+      const IntegrationRule &ir = *irs[mesh.GetElementType(i)];
+      for (int j = 0; j < ir.GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir.IntPoint(j);
+         tr->SetIntPoint(&ip);
+         double val = coeff.Eval(*tr, ip);
+
+         norm += ip.weight * tr->Weight() * val;
+      }
+   }
+   return norm;
+}
+
+#ifdef MFEM_USE_MPI
+double ComputeGlobalMean(Coefficient &coeff, ParMesh &pmesh,
+                         const IntegrationRule *irs[])
+{
+   double loc_norm = ComputeMean(coeff, pmesh, irs);
+   double glob_norm = 0;
+
+   MPI_Comm comm = pmesh.GetComm();
+
+   MPI_Allreduce(&loc_norm, &glob_norm, 1, MPI_DOUBLE, MPI_SUM, comm);
+
+   return glob_norm;
+}
+#endif
+
 }

fem/coefficient.hpp

@@ -2544,6 +2544,19 @@ real_t ComputeGlobalLpNorm(real_t p, VectorCoefficient &coeff, ParMesh &pmesh,
                            const IntegrationRule *irs[]);
 #endif
 
+/* HDG */
+/** Compute the mean of a function f.
+    $ | f |_{mean} = ( \int_\Omega f d\Omega) $ */
+double ComputeMean(Coefficient &coeff, Mesh &mesh,
+                   const IntegrationRule *irs[]);
+
+#ifdef MFEM_USE_MPI
+/** Compute the mean of a function f.
+    $ | f |_{mean} = ( \int_\Omega f d\Omega) $ */
+double ComputeGlobalMean(Coefficient &coeff, ParMesh &pmesh,
+                         const IntegrationRule *irs[]);
+#endif
+
 }
 
 #endif

fem/fe/fe_base.cpp

@@ -133,6 +133,20 @@ void FiniteElement::GetTransferMatrix(const FiniteElement &fe,
    MFEM_ABORT("method is not overloaded");
 }
 
+/* HDG */
+void FiniteElement::Project (
+   Coefficient &coeff, FaceElementTransformations &Trans, Vector &dofs) const
+{
+   mfem_error ("FiniteElement::Project (...) (skeleton) is not overloaded !");
+}
+
+/* HDG */
+void FiniteElement::Project (
+   VectorCoefficient &coeff, FaceElementTransformations &Trans, Vector &dofs) const
+{
+   mfem_error ("FiniteElement::Project (...) (skeleton - VectorCoefficient) is not overloaded !");
+}
+
 void FiniteElement::Project(
    Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
 {
@@ -803,6 +817,52 @@ void NodalFiniteElement::GetLocalRestriction(ElementTransformation &Trans,
    R.Threshold(1e-12);
 }
 
+/* HDG */
+void NodalFiniteElement::Project (
+   Coefficient &coeff, FaceElementTransformations &Trans, Vector &dofs) const
+{
+   for (int i = 0; i < dof; i++)
+   {
+      const IntegrationPoint &ip = Nodes.IntPoint(i);
+      // some coefficients expect that Trans.IntPoint is the same
+      // as the second argument of Eval
+      Trans.Face->SetIntPoint(&ip);
+      dofs(i) = coeff.Eval (*Trans.Face, ip);
+      if (map_type == INTEGRAL)
+      {
+         dofs(i) *= Trans.Face->Weight();
+      }
+   }
+}
+
+/* HDG */
+void NodalFiniteElement::Project (
+   VectorCoefficient &vc, FaceElementTransformations &Trans, Vector &dofs) const
+{
+   MFEM_ASSERT(dofs.Size() == vc.GetVDim()*dof, "");
+
+   Vector x (vc.GetVDim());
+
+   for (int i = 0; i < dof; i++)
+   {
+      const IntegrationPoint &ip = Nodes.IntPoint(i);
+      // some coefficients expect that Trans.IntPoint is the same
+      // as the second argument of Eval
+      Trans.Face->SetIntPoint(&ip);
+      vc.Eval (x, *Trans.Face, ip);
+      if (map_type == INTEGRAL)
+      {
+         x *= Trans.Face->Weight();
+      }
+      for (int j = 0; j < x.Size(); j++)
+      {
+         dofs(dof*j+i) = x(j);
+      }
+   }
+}
+
+
+
 void NodalFiniteElement::Project(
    Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
 {

fem/fe/fe_base.hpp

@@ -238,6 +238,8 @@ public:
 };
 
 class ElementTransformation;
+/* HDG */
+class FaceElementTransformations;
 class Coefficient;
 class VectorCoefficient;
 class MatrixCoefficient;
@@ -523,6 +525,15 @@ public:
                                   ElementTransformation &Trans,
                                   DenseMatrix &I) const;
 
+   /* HDG */
+   /** Given a coefficient and a transformation, compute its projection
+       (approximation) in the local finite dimensional skeleton space in terms
+       of the degrees of freedom. */
+   virtual void Project (Coefficient &coeff,
+                         FaceElementTransformations &Trans, Vector &dofs) const;
+   /* HDG */
+   virtual void Project (VectorCoefficient &coeff,
+                         FaceElementTransformations &Trans, Vector &dofs) const;
    /** @brief Given a coefficient and a transformation, compute its projection
        (approximation) in the local finite dimensional space in terms
        of the degrees of freedom. */
@@ -764,6 +775,13 @@ public:
    void Project(VectorCoefficient &vc,
                 ElementTransformation &Trans, Vector &dofs) const override;
 
+   /* HDG */
+   virtual void Project (Coefficient &coeff,
+                         FaceElementTransformations &Trans, Vector &dofs) const override;
+   /* HDG */
+   virtual void Project (VectorCoefficient &vc,
+                         FaceElementTransformations &Trans, Vector &dofs) const override;
+
    // (mc.height x mc.width) @ DOFs -> (Dof x mc.width x mc.height) in dofs
    void ProjectMatrixCoefficient(
       MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const override;

fem/gridfunc.cpp

@@ -446,6 +446,52 @@ void GridFunction::GetNodalValues(int i, Array<real_t> &nval, int vdim) const
    }
 }
 
+/* HDG */
+double GridFunction::GetValueFacet(FaceElementTransformations &T,
+                                   const IntegrationPoint &ip, int vdim, Vector *tr)
+const
+{
+   if (tr)
+   {
+      T.SetIntPoint(&ip);
+      T.Transform(ip, *tr);
+   }
+
+   const FiniteElement * fe = NULL;
+   Array<int> dofs;
+
+   switch (T.ElementType)
+   {
+      case ElementTransformation::FACE:
+      {
+         fe = fes->GetFaceElement(T.ElementNo);
+         fes->GetFaceDofs(T.ElementNo, dofs);
+         break;
+      }
+      case ElementTransformation::BDR_FACE:
+      {
+         fe = fes->GetBE(T.ElementNo);
+         fes->GetBdrElementDofs(T.ElementNo, dofs);
+         break;
+      }
+   }
+
+   fes->DofsToVDofs(vdim-1, dofs);
+   Vector DofVal(dofs.Size()), LocVec;
+   if (fe->GetMapType() == FiniteElement::VALUE)
+   {
+      fe->CalcShape(ip, DofVal);
+   }
+   else
+   {
+      fe->CalcPhysShape(T, DofVal);
+   }
+   GetSubVector(dofs, LocVec);
+
+   return (DofVal * LocVec);
+
+}
+
 real_t GridFunction::GetValue(int i, const IntegrationPoint &ip, int vdim)
 const
 {
@@ -2261,6 +2307,42 @@ void GridFunction::AccumulateAndCountBdrTangentValues(
    }
 }
 
+/* HDG */
+double GridFunction::ComputeMean(const IntegrationRule *irs[]) const
+{
+   double global_mean = 0.0;
+   double element_mean;
+   const FiniteElement *fe;
+   ElementTransformation *T;
+   Vector vals;
+
+   for (int i = 0; i < fes->GetNE(); i++)
+   {
+      fe = fes->GetFE(i);
+      const IntegrationRule *ir;
+      if (irs)
+      {
+         ir = irs[fe->GetGeomType()];
+      }
+      else
+      {
+         int intorder = 2*fe->GetOrder() + 1; // <----------
+         ir = &(IntRules.Get(fe->GetGeomType(), intorder));
+      }
+      GetValues(i, *ir, vals);
+      T = fes->GetElementTransformation(i);
+      element_mean = 0.0;
+      for (int j = 0; j < ir->GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(j);
+         T->SetIntPoint(&ip);
+         element_mean += ip.weight * T->Weight() * vals(j);
+      }
+      global_mean += element_mean;
+   }
+   return global_mean;
+}
+
 void GridFunction::ComputeMeans(AvgType type, Array<int> &zones_per_vdof)
 {
    switch (type)
@@ -2352,6 +2434,133 @@ void GridFunction::ProjectDeltaCoefficient(DeltaCoefficient &delta_coeff,
    }
 }
 
+/* HDG */
+void GridFunction::ProjectCoefficientSkeleton(Coefficient &coeff)
+{
+   DeltaCoefficient *delta_c = dynamic_cast<DeltaCoefficient *>(&coeff);
+
+   // check if the GridFunction belongs to a skeletal FEColl
+   MFEM_ASSERT(bool(dynamic_cast<const DG_Interface_FECollection*>
+                    (fes->FEColl())) || bool (dynamic_cast<const H1_Trace_FECollection*>
+                                              (fes->FEColl())),
+               "Incorrect FEColl");
+
+   if (delta_c == NULL)
+   {
+      Mesh *mesh = fes->GetMesh();
+      int nfaces = mesh->GetNumFaces();
+
+      Array<int> vdofs;
+      Vector vals;
+
+      for (int i = 0; i < nfaces; i++)
+      {
+         fes->GetFaceVDofs(i, vdofs);
+         vals.SetSize(vdofs.Size());
+         fes->GetFaceElement(i)->Project(coeff, *mesh->GetFaceElementTransformations(i),
+                                         vals);
+         SetSubVector(vdofs, vals);
+      }
+   }
+   else
+   {
+      real_t integral;
+
+      ProjectDeltaCoefficient(*delta_c, integral);
+
+      (*this) *= (delta_c->Scale() / integral);
+   }
+}
+
+void GridFunction::ProjectCoefficientSkeleton(VectorCoefficient &vcoeff)
+{
+   Mesh *mesh = fes->GetMesh();
+   int nfaces = mesh->GetNumFaces();
+   Array<int> vdofs;
+   Vector vals;
+
+   // check if the GridFunction belongs to a skeletal FEColl
+   MFEM_ASSERT(bool(dynamic_cast<const DG_Interface_FECollection*>
+                    (fes->FEColl())) || bool (dynamic_cast<const H1_Trace_FECollection*>
+                                              (fes->FEColl())),
+               "Incorrect FEColl");
+
+   for (int i = 0; i < nfaces; i++)
+   {
+      fes->GetFaceVDofs(i, vdofs);
+      vals.SetSize(vdofs.Size());
+      fes->GetFaceElement(i)->Project(vcoeff, *mesh->GetFaceElementTransformations(i),
+                                      vals);
+      SetSubVector(vdofs, vals);
+   }
+}
+
+void GridFunction::ProjectCoefficientSkeletonBdr(Coefficient &coeff)
+{
+   DeltaCoefficient *delta_c = dynamic_cast<DeltaCoefficient *>(&coeff);
+
+   // check if the GridFunction belongs to a skeletal FEColl
+   MFEM_ASSERT(bool(dynamic_cast<const DG_Interface_FECollection*>
+                    (fes->FEColl())) || bool (dynamic_cast<const H1_Trace_FECollection*>
+                                              (fes->FEColl())),
+               "Incorrect FEColl");
+
+   if (delta_c == NULL)
+   {
+      Mesh *mesh = fes->GetMesh();
+      int nbdrfaces = mesh->GetNBE();
+      Array<int> vdofs;
+      Vector vals;
+
+      for (int i = 0; i < nbdrfaces; i++)
+      {
+         int face = mesh->GetBdrElementFaceIndex(i);
+         fes->GetFaceVDofs(face, vdofs);
+         vals.SetSize(vdofs.Size());
+         fes->GetFaceElement(face)->Project(coeff,
+                                            *mesh->GetFaceElementTransformations(face),
+                                            vals);
+         SetSubVector(vdofs, vals);
+      }
+   }
+   else
+   {
+      real_t integral;
+
+      ProjectDeltaCoefficient(*delta_c, integral);
+
+      (*this) *= (delta_c->Scale() / integral);
+   }
+}
+
+void GridFunction::ProjectCoefficientSkeletonBdr(VectorCoefficient &vcoeff)
+{
+   Mesh *mesh = fes->GetMesh();
+   int nbdrfaces = mesh->GetNBE();
+   Array<int> vdofs;
+   Vector vals;
+
+   // check if the GridFunction belongs to a skeletal FEColl
+   MFEM_ASSERT(bool(dynamic_cast<const DG_Interface_FECollection*>
+                    (fes->FEColl())) || bool (dynamic_cast<const H1_Trace_FECollection*>
+                                              (fes->FEColl())),
+               "Incorrect FEColl");
+
+   for (int i = 0; i < nbdrfaces; i++)
+   {
+      int face = mesh->GetBdrElementFaceIndex(i);
+      fes->GetFaceVDofs(face, vdofs);
+      vals.SetSize(vdofs.Size());
+      fes->GetFaceElement(face)->Project(vcoeff,
+                                         *mesh->GetFaceElementTransformations(face),
+                                         vals);
+      SetSubVector(vdofs, vals);
+   }
+
+}
+
+/* HDG ends */
+
 void GridFunction::ProjectCoefficient(Coefficient &coeff)
 {
    DeltaCoefficient *delta_c = dynamic_cast<DeltaCoefficient *>(&coeff);
@@ -3365,6 +3574,69 @@ real_t GridFunction::ComputeW11Error(
    return error;
 }
 
+/// To compute \| mean(u) - mean(u_h) \|_p
+double GridFunction::ComputeMeanLpError(const double p, Coefficient &exsol,
+                                        const IntegrationRule *irs[]) const
+{
+   double error = 0.0;
+   double aux, local_error, local_size ;
+   const FiniteElement *fe;
+   ElementTransformation *T;
+   Vector vals;
+
+   for (int i = 0; i < fes->GetNE(); i++)
+   {
+      fe = fes->GetFE(i);
+      const IntegrationRule *ir;
+      if (irs)
+      {
+         ir = irs[fe->GetGeomType()];
+      }
+      else
+      {
+         int intorder = 2*fe->GetOrder() + 1; // <----------
+         ir = &(IntRules.Get(fe->GetGeomType(), intorder));
+      }
+      GetValues(i, *ir, vals);
+      T = fes->GetElementTransformation(i);
+      local_error = local_size = 0.0 ;
+      for (int j = 0; j < ir->GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(j);
+         T->SetIntPoint(&ip);
+         aux = (vals(j) - exsol.Eval(*T, ip));
+         local_error += ip.weight * T->Weight() * aux;
+         local_size += ip.weight * T->Weight();
+      }
+
+      if (p < numeric_limits<double>::infinity())
+      {
+         aux = pow(fabs(local_error), p) / pow(local_size, p-1.);
+         error += aux;
+      }
+      else
+      {
+         aux = fabs(local_error) / fabs(local_size);
+         error = std::max(error, aux);
+      }
+   }
+
+   if (p < numeric_limits<double>::infinity())
+   {
+      // negative quadrature weights may cause the error to be negative
+      if (error < 0.)
+      {
+         error = -pow(-error, 1./p);
+      }
+      else
+      {
+         error = pow(error, 1./p);
+      }
+   }
+
+   return error;
+}
+
 real_t GridFunction::ComputeLpError(const real_t p, Coefficient &exsol,
                                     Coefficient *weight,
                                     const IntegrationRule *irs[],

fem/gridfunc.hpp

@@ -23,6 +23,8 @@
 #include <limits>
 #include <ostream>
 #include <string>
+/* HDG */
+//#include "lininteg.hpp"
 
 namespace mfem
 {
@@ -286,6 +288,12 @@ public:
                            DenseMatrix &vals, DenseMatrix &tr) const;
    ///@}
 
+   /* HDG */
+   double GetValueFacet(FaceElementTransformations &T,
+                        const IntegrationPoint &ip,
+                        int vdim = 1,
+                        Vector *tr = NULL) const;
+
    void GetLaplacians(int i, const IntegrationRule &ir, Vector &laps,
                       int vdim = 1) const;
 
@@ -417,6 +425,14 @@ public:
        projection matrix. */
    void ProjectGridFunction(const GridFunction &src);
 
+   /* HDG */
+   void ProjectCoefficientSkeleton(Coefficient &coeff);
+   void ProjectCoefficientSkeleton(VectorCoefficient &vcoeff);
+   //   void ProjectCoefficientSkeletonBdr(Coefficient &coeff,
+   //       Array<int> &bdr_attr_marker);
+   void ProjectCoefficientSkeletonBdr(Coefficient &coeff);
+   void ProjectCoefficientSkeletonBdr(VectorCoefficient &vcoeff);
+
    /** @brief Project @a coeff Coefficient to @a this GridFunction. The
        projection computation depends on the choice of the FiniteElementSpace
        #fes. Note that this is usually interpolation at the degrees of freedom
@@ -493,6 +509,10 @@ protected:
                                            const Array<int> &bdr_attr,
                                            Array<int> &values_counter);
 
+   /* HDG */
+   /* Compute the mean of a GridFunction */
+   double ComputeMean(const IntegrationRule *irs[]) const;
+
    // Complete the computation of averages; called e.g. after
    // AccumulateAndCountZones().
    void ComputeMeans(AvgType type, Array<int> &zones_per_vdof);
@@ -671,6 +691,7 @@ public:
    ///       function uses the absolute values of the element-wise integrals.
    ///       This may lead to results which are not entirely consistent with
    ///       such integration rules.
+
    virtual real_t ComputeGradError(VectorCoefficient *exgrad,
                                    const IntegrationRule *irs[] = NULL) const;
 
@@ -1137,6 +1158,7 @@ public:
    ///       function uses the absolute values of the element-wise integrals.
    ///       This may lead to results which are not entirely consistent with
    ///       such integration rules.
+
    virtual real_t ComputeLpError(const real_t p, Coefficient &exsol,
                                  Coefficient *weight = NULL,
                                  const IntegrationRule *irs[] = NULL,
@@ -1492,6 +1514,12 @@ public:
                                        ) const
    { ComputeElementLpErrors(infinity(), exsol, error, NULL, NULL, irs); }
 
+   /* HDG */
+   double ComputeMeanLpError(const double p, Coefficient &exsol,
+                             const IntegrationRule *irs[] = NULL) const;
+
+
+
    virtual void ComputeFlux(BilinearFormIntegrator &blfi,
                             GridFunction &flux,
                             bool wcoef = true, int subdomain = -1);

fem/linearform.cpp

@@ -105,6 +105,25 @@ void LinearForm::AddInteriorFaceIntegrator(LinearFormIntegrator *lfi)
    interior_face_integs.Append(lfi);
 }
 
+/* HDG */
+void LinearForm::AddSktBoundaryNeumannIntegrator(LinearFormIntegrator * lfi)
+{
+   bdrsklneufi.Append (lfi);
+   bdrsklneufi_marker.Append(NULL); // NULL -> all attributes are active
+}
+void LinearForm::AddSktBoundaryNeumannIntegrator(LinearFormIntegrator * lfi,
+                                                 Array<int> &bdr_attr_marker)
+{
+   bdrsklneufi.Append (lfi);
+   bdrsklneufi_marker.Append(
+      &bdr_attr_marker); // NULL -> all attributes are active
+}
+/* HDG */
+void LinearForm::AddSktInteriorFaceIntegrator(LinearFormIntegrator * lfi)
+{
+   interiorsklfi.Append (lfi);
+}
+
 bool LinearForm::SupportsDevice() const
 {
    // return false for NURBS meshes, so we don’t convert it to non-NURBS
@@ -338,6 +357,77 @@ void LinearForm::Assemble()
          }
       }
    }
+
+   /* HDG */
+   if (bdrsklneufi.Size())
+   {
+      FaceElementTransformations *ftr;
+      const FiniteElement *face_fe;
+      int nbdrfaces = fes->GetNBE();
+      Mesh *mesh = fes -> GetMesh();
+
+      for (int i = 0; i < nbdrfaces; i++)
+      {
+         int face = mesh->GetBdrElementFaceIndex(i);
+         ftr = mesh->GetBdrFaceTransformations(
+                  i); // the transformation of the face
+         //          fes->GetBdrElementVDofs(i, vdofs);   // the degrees of freedom related to the face
+         fes->GetFaceVDofs(face, vdofs);   // the degrees of freedom related to the face
+         face_fe = fes->GetFaceElement(
+                      face);   // point face_fe to the FiniteElement over the edge
+
+         if (ftr != NULL)
+         {
+            for (int k = 0; k < bdrsklneufi.Size(); k++) // Loop over the related interals
+            {
+               int compute = 0;
+               if (bdrsklneufi_marker[k] == NULL)
+               {
+                  compute = 1;
+               }
+               else
+               {
+                  Array<int> &bdr_marker = *bdrsklneufi_marker[k];
+                  const int bdr_attr = mesh->GetBdrAttribute(i);
+                  if (bdr_marker[bdr_attr-1] == 1)
+                  {
+                     compute = 1;
+                  }
+               }
+               if (compute)
+               {
+                  bdrsklneufi[k] -> AssembleRHSElementVect (*face_fe, *ftr, elemvect);
+                  AddElementVector (vdofs, elemvect);
+               }
+            }
+         }
+
+      }
+   }
+
+   /* HDG */
+   if (interiorsklfi.Size())
+   {
+      Mesh *mesh = fes->GetMesh();
+
+      for (int k = 0; k < interiorsklfi.Size(); k++)
+      {
+         for (int i = 0; i < mesh->GetNumFaces(); i++)
+         {
+            FaceElementTransformations *tr = NULL;
+            tr = mesh->GetInteriorFaceTransformations (i);
+            if (tr != NULL)
+            {
+               fes->GetFaceVDofs(i, vdofs);   // the degrees of freedom related to the face
+
+               interiorsklfi[k]->
+               AssembleRHSElementVect(*fes->GetFaceElement(i),
+                                      *tr, elemvect);
+               AddElementVector (vdofs, elemvect);
+            }
+         }
+      }
+   }
 }
 
 void LinearForm::Update()
@@ -428,6 +518,8 @@ LinearForm::~LinearForm()
       { delete boundary_face_integs[k]; }
       for (k=0; k < interior_face_integs.Size(); k++)
       { delete interior_face_integs[k]; }
+      for (k=0; k < bdrsklneufi.Size(); k++)
+      { delete bdrsklneufi[k]; }
    }
 
    delete ext;

fem/linearform.hpp

@@ -71,6 +71,16 @@ protected:
    /// The reference coordinates where the centers of the delta functions lie
    Array<IntegrationPoint> domain_delta_integs_ip;
 
+   /* HDG */
+   /// Set of Boundary Skeleton Integrators .
+   Array<LinearFormIntegrator*> bdrsklneufi;
+   Array<Array<int>*>           bdrsklneufi_marker;
+
+   /* HDG */
+   /// Set of Interior Skeleton Integrators .
+   Array<LinearFormIntegrator*> interiorsklfi;
+
+
    /// If true, the delta locations are not (re)computed during assembly.
    bool HaveDeltaLocations()
    { return (domain_delta_integs_elem_id.Size() != 0); }
@@ -151,6 +161,17 @@ public:
    /// Adds new Boundary Face Integrator. Assumes ownership of @a lfi.
    void AddBdrFaceIntegrator(LinearFormIntegrator *lfi);
 
+   /* HDG */
+   /// Adds new Boundary Face Integrator with face number.
+   void AddSktBoundaryNeumannIntegrator (LinearFormIntegrator * lfi);
+
+   void AddSktBoundaryNeumannIntegrator (LinearFormIntegrator * lfi,
+                                         Array<int> &bdr_attr_marker);
+
+   /* HDG */
+   /// Access all integrators added with AddSktBoundaryNeumannIntegrator().
+   Array<LinearFormIntegrator*> *GetBDRSKTFLFI() { return &bdrsklneufi; }
+
    /** @brief Add new Boundary Face Integrator, restricted to the given boundary
        attributes.
 
@@ -190,6 +211,13 @@ public:
        corresponding pointer (to Array<int>) will be NULL. */
    Array<Array<int>*> *GetFLFI_Marker() { return &boundary_face_integs_marker; }
 
+   /* HDG */
+   void AddSktInteriorFaceIntegrator(LinearFormIntegrator * lfi);
+
+   /* HDG */
+   /// Access all integrators added with AddSktInteriorFaceIntegrator().
+   Array<LinearFormIntegrator*> *GetISKTFLFI() { return &interiorsklfi; }
+
    /// @brief Which assembly algorithm to use: the new device-compatible fast
    /// assembly (true), or the legacy CPU-only algorithm (false).
    /** If not set, the default value is false.  If used, this method must be

fem/lininteg.cpp

@@ -22,6 +22,13 @@ void LinearFormIntegrator::AssembleDevice(const FiniteElementSpace &fes,
    MFEM_ABORT("Not supported.");
 }
 
+
+void LinearFormIntegrator::AssembleRHSElementVect(
+   const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
+{
+   mfem_error("LinearFormIntegrator::AssembleRHSElementVect(..., ElementTransformation, ...)");
+}
+
 void LinearFormIntegrator::AssembleRHSElementVect(
    const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
 {
@@ -121,6 +128,8 @@ void DomainLFGradIntegrator::AssembleRHSElementVect(
    }
 }
 
+
+
 void DomainLFGradIntegrator::AssembleDeltaElementVect(
    const FiniteElement &fe, ElementTransformation &Trans, Vector &elvect)
 {

fem/lininteg.hpp

@@ -56,7 +56,6 @@ public:
    virtual ~LinearFormIntegrator() { }
 };
 
-
 /// Abstract class for integrators that support delta coefficients
 class DeltaLFIntegrator : public LinearFormIntegrator
 {

fem/pbilinearform.cpp

@@ -262,9 +262,39 @@ void ParBilinearForm::AssembleSharedFaces(int skip_zeros)
    }
 }
 
+/* HDG */
+void ParBilinearForm::AssembleSharedHDGFaces(int skip_zeros)
+{
+   ParMesh *pmesh = pfes->GetParMesh();
+   FaceElementTransformations *ftr;
+   const FiniteElement *face_fe;
+   Array<int> vdofs;
+   DenseMatrix elemmat;
+
+   int nfaces = pmesh->GetNSharedFaces();
+   for (int i = 0; i < nfaces; i++)
+   {
+      int FaceNo = pmesh->GetSharedFace(i);
+      ftr = pmesh->GetSharedFaceTransformations(i);
+      pfes->GetFaceVDofs(FaceNo, vdofs);
+      face_fe = pfes->GetFaceElement(FaceNo);
+      //face_fe = pfes->GetFaceNbrFaceFE(FaceNo);
+      if (ftr != NULL)
+      {
+         for (int k = 0; k < hdgintbfi.Size(); k++)
+         {
+            hdgintbfi[k]->AssembleFaceMatrix(*face_fe, *ftr, elemmat);
+            elemmat *= 0.5;
+            mat->AddSubMatrix(vdofs, vdofs, elemmat, skip_zeros);
+         }
+      }
+   }
+}
+
 void ParBilinearForm::Assemble(int skip_zeros)
 {
-   if (interior_face_integs.Size())
+   /* HDG */
+   if (interior_face_integs.Size() || hdgintbfi.Size())
    {
       pfes->ExchangeFaceNbrData();
       if (!ext && mat == NULL)
@@ -279,6 +309,13 @@ void ParBilinearForm::Assemble(int skip_zeros)
    {
       AssembleSharedFaces(skip_zeros);
    }
+
+   /* HDG */
+   if (!ext && hdgintbfi.Size() > 0)
+   {
+      AssembleSharedHDGFaces(skip_zeros);
+   }
+
 }
 
 void ParBilinearForm::AssembleDiagonal(Vector &diag) const

fem/pbilinearform.hpp

@@ -43,6 +43,8 @@ protected:
 
    void AssembleSharedFaces(int skip_zeros = 1);
 
+   /* HDG */
+   void AssembleSharedHDGFaces(int skip_zeros = 1);
 private:
    /// Copy construction is not supported; body is undefined.
    ParBilinearForm(const ParBilinearForm &);

fem/pgridfunc.cpp

@@ -565,6 +565,50 @@ void ParGridFunction::ProjectCoefficient(Coefficient &coeff)
    }
 }
 
+/* HDG */
+void ParGridFunction::ProjectCoefficientSkeleton(Coefficient &coeff)
+{
+   DeltaCoefficient *delta_c = dynamic_cast<DeltaCoefficient *>(&coeff);
+
+   if (delta_c == NULL)
+   {
+      GridFunction::ProjectCoefficientSkeleton(coeff);
+   }
+   else
+   {
+      real_t loc_integral, glob_integral;
+
+      ProjectDeltaCoefficient(*delta_c, loc_integral);
+
+      MPI_Allreduce(&loc_integral, &glob_integral, 1, MPI_DOUBLE, MPI_SUM,
+                    pfes->GetComm());
+
+      (*this) *= (delta_c->Scale() / glob_integral);
+   }
+}
+
+/* HDG */
+void ParGridFunction::ProjectCoefficientSkeletonBdr(Coefficient &coeff)
+{
+   DeltaCoefficient *delta_c = dynamic_cast<DeltaCoefficient *>(&coeff);
+
+   if (delta_c == NULL)
+   {
+      GridFunction::ProjectCoefficientSkeletonBdr(coeff);
+   }
+   else
+   {
+      real_t loc_integral, glob_integral;
+
+      ProjectDeltaCoefficient(*delta_c, loc_integral);
+
+      MPI_Allreduce(&loc_integral, &glob_integral, 1, MPI_DOUBLE, MPI_SUM,
+                    pfes->GetComm());
+
+      (*this) *= (delta_c->Scale() / glob_integral);
+   }
+}
+
 void ParGridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff)
 {
    // local maximal element attribute for each dof
@@ -1215,6 +1259,15 @@ void ParGridFunction::SaveAsOne(std::ostream &os) const
    delete [] nrdofs;
 }
 
+/* HDG */
+double GlobalMean(double loc_mean, MPI_Comm comm)
+{
+   double glob_mean;
+   MPI_Allreduce(&loc_mean, &glob_mean, 1, MPI_DOUBLE, MPI_SUM, comm);
+
+   return glob_mean;
+}
+
 real_t GlobalLpNorm(const real_t p, real_t loc_norm, MPI_Comm comm)
 {
    real_t glob_norm;

fem/pgridfunc.hpp

@@ -45,6 +45,11 @@ namespace mfem
 ///       such integration rules.
 real_t GlobalLpNorm(const real_t p, real_t loc_norm, MPI_Comm comm);
 
+/* HDG */
+/* Compute the mean of a coefficient in parallel */
+double GlobalMean(double loc_mean, MPI_Comm comm);
+
+
 /// Class for parallel grid function
 class ParGridFunction : public GridFunction
 {
@@ -259,6 +264,12 @@ public:
    using GridFunction::ProjectCoefficient;
    void ProjectCoefficient(Coefficient &coeff) override;
 
+   /* HDG */
+   using GridFunction::ProjectCoefficientSkeleton;
+   void ProjectCoefficientSkeleton(Coefficient &coeff);
+   using GridFunction::ProjectCoefficientSkeletonBdr;
+   void ProjectCoefficientSkeletonBdr(Coefficient &coeff);
+
    using GridFunction::ProjectDiscCoefficient;
    /** @brief Project a discontinuous vector coefficient as a grid function on
        a continuous finite element space. The values in shared dofs are
@@ -433,6 +444,7 @@ public:
    ///       function uses the absolute values of the element-wise integrals.
    ///       This may lead to results which are not entirely consistent with
    ///       such integration rules.
+
    real_t ComputeDGFaceJumpError(Coefficient *exsol,
                                  Coefficient *ell_coeff,
                                  JumpScaling jump_scaling,
@@ -577,6 +589,12 @@ public:
                              p, exsol, weight, v_weight, irs), pfes->GetComm());
    }
 
+   /* HDG */
+   double ComputeMean(const IntegrationRule *irs[] = NULL) const
+   {
+      return GlobalMean(GridFunction::ComputeMean(irs), pfes->GetComm());
+   }
+
    void ComputeFlux(BilinearFormIntegrator &blfi,
                     GridFunction &flux,
                     bool wcoef = true, int subdomain = -1) override;

fem/plinearform.cpp

@@ -52,6 +52,11 @@ void ParLinearForm::Assemble()
       pfes->ExchangeFaceNbrData();
       AssembleSharedFaces();
    }
+   if (interiorsklfi.Size())
+   {
+      pfes->ExchangeFaceNbrData();
+      AssembleSharedHDGFaces();
+   }
 }
 
 bool ParLinearForm::SupportsDevice() const
@@ -93,6 +98,37 @@ void ParLinearForm::AssembleSharedFaces()
    }
 }
 
+void ParLinearForm::AssembleSharedHDGFaces()
+{
+   Array<int> vdofs;
+   Vector elemvect;
+
+   if (interiorsklfi.Size())
+   {
+      ParMesh *pmesh = pfes->GetParMesh();
+
+      for (int k = 0; k < interiorsklfi.Size(); k++)
+      {
+         for (int i = 0; i < pmesh->GetNSharedFaces(); i++)
+         {
+            FaceElementTransformations *tr = NULL;
+            tr = pmesh->GetSharedFaceTransformations(i);
+            if (tr != NULL)
+            {
+               int face_idx = pmesh->GetSharedFace(i);
+               fes->GetFaceVDofs(face_idx,
+                                 vdofs);   // the degrees of freedom related to the face
+
+               interiorsklfi[k]->
+               AssembleRHSElementVect(*fes->GetFaceElement(face_idx),
+                                      *tr, elemvect);
+               AddElementVector (vdofs, elemvect);
+            }
+         }
+      }
+   }
+}
+
 void ParLinearForm::ParallelAssemble(Vector &tv)
 {
    const Operator* prolong = pfes->GetProlongationMatrix();

fem/plinearform.hpp

@@ -124,6 +124,8 @@ public:
 
    void AssembleSharedFaces();
 
+   void AssembleSharedHDGFaces();
+
    /// Assemble the vector on the true dofs, i.e. P^t v.
    void ParallelAssemble(Vector &tv);
 

makefile

@@ -125,7 +125,8 @@ EXAMPLE_TEST_DIRS := examples
 
 MINIAPP_SUBDIRS = common electromagnetics meshing navier performance tools \
  toys nurbs gslib adjoint solvers shifted mtop parelag tribol autodiff dfem \
- hooke multidomain dpg hdiv-linear-solver spde diag-smoothers
+ hooke multidomain dpg hdiv-linear-solver spde diag-smoothers hdg
+
 MINIAPP_DIRS := $(addprefix miniapps/,$(MINIAPP_SUBDIRS))
 MINIAPP_TEST_DIRS := $(filter-out %/common,$(MINIAPP_DIRS))
 MINIAPP_USE_COMMON := $(addprefix miniapps/,electromagnetics meshing tools \

mesh/mesh.cpp

@@ -1616,7 +1616,7 @@ Element::Type Mesh::GetFaceElementType(int Face) const
 
 Array<int> Mesh::GetFaceToBdrElMap() const
 {
-   Array<int> face_to_be(Dim == 2 ? NumOfEdges : NumOfFaces);
+   Array<int> face_to_be(GetNumFaces());
    face_to_be = -1;
    for (int i = 0; i < NumOfBdrElements; i++)
    {

mesh/mesh.hpp

@@ -1604,6 +1604,10 @@ public:
    void GetElementVertices(int i, Array<int> &v) const
    { elements[i]->GetVertices(v); }
 
+   /// HDG:sets the indices of the vertices of element i.
+   void SetElementVertices(int i, Array<int> &v) const
+   { elements[i]->SetVertices(v); }
+
    /// Returns the indices of the vertices of boundary element i.
    void GetBdrElementVertices(int i, Array<int> &v) const
    { boundary[i]->GetVertices(v); }
@@ -1936,7 +1940,6 @@ public:
                                   IsoparametricTransformation &ElTr2) const;
 
    /// @}
-
    /// @anchor mfem_Mesh_geom_factors
    /// @name Access the coordinate transformation at quadrature points
    ///

miniapps/hdg/HDGBilinearForm.hpp

@@ -0,0 +1,330 @@
+// 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.
+//
+// Implementation of class HDGBilinearForm
+//
+// Contributed by: T. Horvath: Oakland University
+//                 S. Rhebergen, A. Sivas: University of Waterloo
+
+#ifndef MFEM_HDGBILINEARFORM
+#define MFEM_HDGBILINEARFORM
+
+#include "../../config/config.hpp"
+#include "../../linalg/linalg.hpp"
+#include "../../fem/fespace.hpp"
+#include "../../fem/gridfunc.hpp"
+#include "../../fem/linearform.hpp"
+#include "../../fem/bilininteg.hpp"
+
+#ifdef MFEM_USE_MPI
+#include <mpi.h>
+#endif
+
+using namespace std;
+using namespace mfem;
+
+class HDGBilinearForm
+{
+protected:
+   /// FE spaces on which the form lives.
+   Array<FiniteElementSpace*> volume_fes, skeletal_fes;
+
+   int NVolumeFES, NSkeletalFES;
+
+   bool parallel;
+
+   /// Sparse matrix to be assembled
+   Array<SparseMatrix*> mat;
+
+   /// Right hand side vector to be assembled.
+   Array<Vector*> rhs_SC;
+
+   /// Table that contains the faces for all elements
+   Table *el_to_face;
+
+   /// List that separates the interior edges from the shared edges
+   Array<int> ess_dofs, Edge_to_SharedEdge;
+
+   /// HDG Integrators
+   Array<BilinearFormIntegrator*> hdg_dbfi;
+   Array<BilinearFormIntegrator*> hdg_fbfi;
+
+   /// Dense matrices to be used for computing the integrals
+   DenseMatrix elemmat1, elemmat2, elemmat3, elemmat4;
+
+   /// Vectors to store A and B, the corresponding offsets and the number
+   /// of elements on which A and B will be stored
+   Array<int> A_offsets, B_offsets;
+   real_t *A_data, *B_data;
+   int elements_A, elements_B;
+
+
+   // may be used in the construction of derived classes
+   HDGBilinearForm()
+   {
+      for (int i =0; i<NVolumeFES; i++)
+      {
+         delete volume_fes[i];
+      }
+      for (int i =0; i<NSkeletalFES; i++)
+      {
+         delete skeletal_fes[i];
+         delete rhs_SC[i];
+      }
+      for (int i =0; i<NSkeletalFES*NSkeletalFES; i++)
+      {
+         delete mat[i];
+      }
+      NVolumeFES = 0;
+      NSkeletalFES = 0;
+      volume_fes = NULL;
+      skeletal_fes = NULL;
+      parallel = false;
+      el_to_face = NULL;
+      A_data = NULL;
+      B_data = NULL;
+      elements_A = elements_B = 0;
+   }
+
+public:
+   /// Creates bilinear form associated with FE spaces *_fes1 and _fes2.
+   HDGBilinearForm(Array<FiniteElementSpace*> &_fes1,
+                   Array<FiniteElementSpace*> &_fes2,
+                   bool _parallel = false);
+
+   // Advection-reaction test case without FES arrays
+   HDGBilinearForm(FiniteElementSpace *_fes1,
+                   FiniteElementSpace *_fes2,
+                   bool _parallel = false);
+
+   // Diffusion test case without FES arrays
+   HDGBilinearForm(FiniteElementSpace *_fes1,
+                   FiniteElementSpace *_fes2,
+                   FiniteElementSpace *_fes3,
+                   bool _parallel = false);
+
+   // Arrays of the HDG domain integrators
+   Array<BilinearFormIntegrator*> *GetHDG_DBFI()
+   {
+      return &hdg_dbfi;
+   }
+
+   // Arrays of the HDG face integrators
+   Array<BilinearFormIntegrator*> *GetHDG_FBFI()
+   {
+      return &hdg_fbfi;
+   }
+
+   /// Finalizes the matrix
+   virtual void Finalize(int skip_zeros = 1);
+
+   // Gets all dofs for a given element (goes over all volume FES)
+   void GetInteriorVDofs(int i, Array<int> &vdofs) const;
+
+   // Gets the values of volume GFs at all dofs for a given element (goes over all volume FES)
+   void GetInteriorSubVector(const Array<GridFunction*> &rhs_gridfunctions,
+                             int i, int ndof, Vector &SubVector) const;
+
+   // Gets the values of a skeletal GFs at all dofs for a given element (goes over all skeletal FES)
+   void GetFaceSubVector(const Array<GridFunction*> &face_gridfunctions,
+                         int i, int ndof, Vector &SubVector) const;
+
+   // Sets the values of volume GFs at all dofs for a given element (goes over all volume FES)
+   void SetInteriorSubVector(Array<GridFunction*> &sol_gridfunctions,
+                             int i, int ndof, Vector &SubVector);
+
+   // Gets all dofs for a given element (goes over all skeletal FES)
+   // includes an array that counts the number of dofs per FES
+   void GetFaceVDofs(int i, Array<int> &vdofs, Array<int> &dof_length) const;
+
+   // Gets all dofs for a given element (goes over all volume FES)
+   void GetFaceVDofs(int i, Array<int> &vdofs) const;
+
+   /// Returns the sparse Schur complement matrix
+   // For block systems: the matrices are counted row-wise
+   const SparseMatrix *SpMatSC(int m = 0) const
+   {
+      MFEM_VERIFY(mat[m], "mat is NULL and can't be dereferenced");
+      return mat[m];
+   }
+
+   /// Returns the
+   // For block systems: the matrices are counted row-wise
+   SparseMatrix *SpMatSC(int m = 0)
+   {
+      MFEM_VERIFY(mat[m], "mat is NULL and can't be dereferenced");
+      return mat[m];
+   }
+
+   /// Returns a constant reference to the right hand side vector
+   const Vector *VectorSC(int m = 0) const
+   {
+      return rhs_SC[m];
+   }
+
+   /// Returns a reference to the right hand side vector
+   Vector *VectorSC(int m = 0)
+   {
+      return rhs_SC[m];
+   }
+
+
+   /// Adds new HDG Integrators (domain terms).
+   void AddHDGDomainIntegrator(BilinearFormIntegrator *bfi);
+
+   /// Adds new HDG Integrators (facet terms).
+   void AddHDGFaceIntegrator(BilinearFormIntegrator *bfi);
+
+   /// Allocates the vectors for the part of A and B that will be stored
+   void Allocate(const Array<int> &bdr_attr_is_ess,
+                 const real_t memA = 0.0, const real_t memB = 0.0);
+
+   /// Assembles the Schur complement - general approach, with FES arrays.
+   // rhs_F is the volume right hand side gf array
+   // rhs_G is the skeleton right hand side gf array
+   // sol is the array of gfs containing the known boundary conditions
+   void AssembleSC(Array<GridFunction*> rhs_F,
+                   Array<GridFunction*> rhs_G,
+                   const Array<int> &bdr_attr_is_ess,
+                   Array<GridFunction*> sol,
+                   int skip_zeros = 1);
+
+   /// Assembles the Schur complement - for the hdg_advection.cpp and hdg_advectionp.cpp test cases
+   void AssembleSC(GridFunction *F,
+                   const real_t memA = 0.0, const real_t memB = 0.0,
+                   int skip_zeros = 1);
+
+   /// Assembles the Schur complement - for the hdg_poisson.cpp and hdg_poissonp.cpp test cases
+   void AssembleSC(GridFunction *F1,
+                   GridFunction *F2,
+                   Array<int> &bdr_attr_is_ess,
+                   GridFunction &sol,
+                   const real_t memA = 0.0, const real_t memB = 0.0,
+                   int skip_zeros = 1);
+
+   /** Eliminate the boundary condition
+    * Currently only used for the diffusion case, but written in a way that it can be used for mutiple skeleton veriables
+    * Eliminates the row and also the column (to keep the SPD matrix for diffusion)
+    */
+   void Eliminate_BC(const Array<int> &vdofs_e1, const Array<int> &vdofs_e1_length,
+                     const int ndof_u,
+                     Array<GridFunction*> sol, Vector *rhs_Volume, Vector *rhs_Skeleton,
+                     DenseMatrix *B_local, DenseMatrix *C_local, DenseMatrix *D_local);
+
+   /* To calculate the inverse of the local matrix A
+    * Can be optimized is some parts are easy to implement (such as the vector mass integrator A11 for the diffusion case)
+    */
+   DenseMatrix CalculateInverse(DenseMatrix A_local);
+
+   /// Computes domain based integrators
+   void compute_domain_integrals(const int elem, DenseMatrix *A_local);
+
+   /// Computes face based integrators
+   void compute_face_integrals(const int elem,
+                               const int edge,
+                               const int isshared,
+                               const bool reconstruct_only,
+                               DenseMatrix *A_local,
+                               DenseMatrix *B_local,
+                               DenseMatrix *C_local,
+                               DenseMatrix *D_local);
+
+   /// Reconstructs u from the facet unknowns - general approach, with FES arrays.
+   // Volume_GF is the volume right hand side gf array
+   // Skeleton_GF is the skeleton solution gf array
+   // u is the volume solution gf array
+   void Reconstruct(Array<GridFunction*> Volume_GF,
+                    Array<GridFunction*> Skeleton_GF,
+                    Array<GridFunction*> u);
+
+   /// Reconstructs u from the facet unknowns - for the hdg_advection.cpp and hdg_advectionp.cpp test cases
+   void Reconstruct(GridFunction *F,
+                    GridFunction *ubar,
+                    GridFunction *u);
+
+   /// Reconstructs u and q from the facet unknowns - for the hdg_poisson.cpp and hdg_poissonp.cpp test cases
+   void Reconstruct(GridFunction *R,
+                    GridFunction *F,
+                    GridFunction *ubar,
+                    GridFunction *q,
+                    GridFunction *u);
+
+   /**
+    Depending on 'assemble' it either assembles the Schur complement, or
+    reconstruct the volume unknowns from the facet solution.
+    For the asssebly part:
+     - Vol_GF is the array of the volume equations right hand side grid functions (might be all zero)
+     - Skel_GF is the array of the skeletal equations right hand side grid functions (might be all zero)
+     - bdr_sol_sol_GF is the projection of the exact solution to the facet unkowns (used only on the boundary for elimination)
+     - bdr_attr_is_ess is the array is essential bouddaries for all facet eqautions (might be empty)
+    For the reconstructions part:
+     - Vol_GF is the array of the volume equations right hand side grid functions (might be all zero)
+     - Skel_GF is the array of the skeletal solutions
+     - bdr_sol_sol_GF is the array of the volume unknowns (to be recontructed)
+
+    */
+   void AssembleReconstruct(Array<GridFunction*> Vol_GF,
+                            Array<GridFunction*> Skel_GF,
+                            const Array<int> &bdr_attr_is_ess,
+                            Array<GridFunction*> bdr_sol_sol_GF,
+                            bool assemble = true,
+                            const real_t memA = 0.0, const real_t memB = 0.0,
+                            int skip_zeros = 1);
+
+
+   /// Updates the spaces
+   virtual void Update(FiniteElementSpace *nfes1 = NULL,
+                       FiniteElementSpace *nfes2 = NULL);
+
+   // Add the vector v_add to the right hand side of the Shur complement system. Works even for block Shur complement systems.
+   void AddToRHS(Array<int> &skeletal_vdofs, Array<int> &skeletal_vdof_length,
+                 Vector v_add);
+
+   // Add the matrix dm_add to the Schur complement system. Works even for block Shur complement systems.
+   void AddToMat(Array<int> &skeletal_vdofs_edge_i,
+                 Array<int> &skeletal_vdof_length_edge_i,
+                 Array<int> &skeletal_vdofs_edge_j, Array<int> &skeletal_vdof_length_edge_j,
+                 DenseMatrix dm_add, int skip_zeros);
+
+   /// Destroys bilinear form.
+   virtual ~HDGBilinearForm();
+
+#ifdef MFEM_USE_MPI
+   // Returns the matrix assembled on the true dofs, i.e. P^t A P.
+   // For block systems: the matrices are counted row-wise
+   HypreParMatrix *ParallelAssembleSC(int i = 0)
+   {
+      return ParallelAssembleSC(i,mat[i]);
+   }
+   // Return the matrix m assembled on the true dofs, i.e. P^t A P
+   HypreParMatrix *ParallelAssembleSC(int i, SparseMatrix *m);
+
+   // Return the right hand side vector
+   HypreParVector *ParallelVectorSC(int i = 0);
+
+   // Return the right hand side vector
+   void ParallelVectorSC(int i, Vector &tv);
+
+   // 2025 Sept begins
+   HypreParMatrix *ParallelAssemble(int i, SparseMatrix *m);
+
+   void ParallelAssemble(OperatorHandle &A, int i = 0) { ParallelAssemble(i, A, mat[i]); }
+
+   /** Returns the matrix @a A_local assembled on the true dofs, i.e.
+      @a A = P^t A_local P in the format (type id) specified by @a A. */
+   // needs input i to know which part of the Schur complement are we working on
+   void ParallelAssemble(int i, OperatorHandle &A, SparseMatrix *m);
+   // 2025 Sept ends
+
+#endif
+};
+
+#endif

miniapps/hdg/hdg_advection.cpp

@@ -0,0 +1,470 @@
+//                        MFEM EDG/HDG example
+//
+// Compile with: make advection
+//
+// Sample runs:  hdg_advection -o 1 -r 1 -tr 4 -no-vis
+//            hdg_advection -o 5 -r 1 -tr 4 -no-vis
+//            hdg_advection -o 1 -r 4 -tr 1
+//            hdg_advection -o 5 -r 4 -tr 1
+//            hdg_advection -o 1 -r 1 -tr 4 -no-vis -m ../data/inline-tri.mesh
+//            hdg_advection -o 5 -r 1 -tr 4 -no-vis -m ../data/inline-tri.mesh
+//            hdg_advection -o 1 -r 5 -tr 1 -m ../data/inline-tri.mesh
+//            hdg_advection -o 5 -r 5 -tr 1 -m ../data/inline-tri.mesh
+//
+//
+// Description:  This example code demonstrates the use of MFEM to define a
+//               finite element discretization of the advection-reaction problem
+//               mu u + a.grad(u) = f with inhomogeneous Neumann boundary conditions.
+//               Specifically, we discretize using a HDG space of the
+//               specified order.
+//
+// The weak form is: seek (u,ubar) such that for all (v, vbar)
+//
+// \mu (u,v)   + (v, a.grad(u) - < 1, [zeta a.n u v] >  + < ubar, [zeta a.n v] > = (f, w)
+// < ubar, [zeta a.n v] >  +
+//       < 1, [zeta a.n ubar vbar] > + < 1, [(1-zeta) a.n ubar vbar >_{\Gamma_N} = < g, vbar >
+//
+// where (.,.) is the d-dimensional L2 product, <.,.> is the d-1 dimensional L2 product,
+// zeta = 1 for inflow boundaries, and 0 otherwise.
+//
+// The discretization is based on the paper:
+//
+// G. N. Wells, Analysis of an interface stabilized finite element method: the advection-diffusion-reaction equation, SIAM J. Numer. Anal., 2011, 49:1, 87--109.
+//
+// Contributed by: T. Horvath, Oakland University
+//                 S. Rhebergen, A. Sivas, University of Waterloo
+
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include <algorithm>
+#include "HDGBilinearForm.hpp"
+#include "hdg_integrators.hpp"
+
+using namespace std;
+using namespace mfem;
+
+//---------------------------------------------------------------------
+// Exact solution and r.h.s.. See below for implementation.
+real_t u_exact(const Vector &x);
+real_t f_rhs  (const Vector &x);
+void advection_function(const Vector &x, Vector &v);
+int dim;
+//---------------------------------------------------------------------
+int main(int argc, char *argv[])
+{
+   StopWatch chrono;
+
+   // 1. Parse command-line options.
+   const char *mesh_file = "../../data/inline-tri.mesh";
+   int order = 1;
+   int total_ref_levels = 2;
+   int initial_ref_levels = 0;
+   bool visualization = true;
+   bool save = true;
+   bool hdg = true;
+   real_t memA = 0.0;
+   real_t memB = 0.0;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree > 1).");
+   args.AddOption(&initial_ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly for the initial calculation.");
+   args.AddOption(&total_ref_levels, "-tr", "--refine",
+                  "Number of times to refine the mesh uniformly.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&save, "-save", "--save-files", "-no-save",
+                  "--no-save-files",
+                  "Enable or disable file saving.");
+   args.AddOption(&hdg, "-hdg", "--hybrid", "-edg",
+                  "--embedded",
+                  "HDG / EDG option.");
+   args.AddOption(&memA, "-memA", "--memoryA",
+                  "Storage of A.");
+   args.AddOption(&memB, "-memB", "--memoryB",
+                  "Storage of B.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   if (order < 1)
+   {
+      cout << "Polynomial order should be > 0. Changing to order 1.";
+      order = 1;
+   }
+
+   // memA, memB \in [0,1], memB <= memA
+   if (memB > memA)
+   {
+      std::cout << "memB cannot be more than memA. Resetting to be equal" << std::endl
+                << std::flush;
+      memA = memB;
+   }
+   if (memA > 1.0)
+   {
+      std::cout << "memA cannot be more than 1. Resetting to 1" << std::endl <<
+                std::flush;
+      memA = 1.0;
+   }
+   else if (memA < 0.0)
+   {
+      std::cout << "memA cannot be less than 0. Resetting to 0." << std::endl <<
+                std::flush;
+      memA = 0.0;
+   }
+   if (memB > 1.0)
+   {
+      std::cout << "memB cannot be more than 1. Resetting to 1" << std::endl <<
+                std::flush;
+      memB = 1.0;
+   }
+   else if (memB < 0.0)
+   {
+      std::cout << "memB cannot be less than 0. Resetting to 0." << std::endl <<
+                std::flush;
+      memB = 0.0;
+   }
+
+   // 2. Read the mesh from the given mesh file. Refine it up to the initial_ref_levels.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+
+   if (mesh->Nonconforming())
+   {
+      cout << "The current implementation does not support Nonconforming meshes. Terminating"
+           << endl << flush;
+      return 1;
+   }
+
+   dim = mesh->Dimension();
+
+   for (int ii=0; ii<initial_ref_levels; ii++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 3. Define the vectors that will contain the errors and the iteration count at every refinement level
+   Vector l2errors(total_ref_levels);
+   Array<int> iterativeMethodIts(total_ref_levels);
+
+   // 4. Define the finite element spaces on the mesh.
+   FiniteElementCollection *Uh_fec(new DG_FECollection(order, dim));
+   FiniteElementCollection *Uhbar_fec = NULL;
+   if (hdg)
+   {
+      Uhbar_fec = new DG_Interface_FECollection(order, dim);
+   }
+   else
+   {
+      Uhbar_fec = new H1_Trace_FECollection(order, dim);
+   }
+
+   // Finite element spaces:
+   // Uh_space is the DG space on elements
+   // ubar_space is the DG space on faces
+   FiniteElementSpace *Uh_space(new FiniteElementSpace(mesh, Uh_fec));
+   FiniteElementSpace *Uhbar_space(new FiniteElementSpace(mesh, Uhbar_fec));
+
+   // 5. Define the coefficients
+   ConstantCoefficient mu(1.0); // reaction constant
+   // Given boundary condition / exact solution
+   FunctionCoefficient ucoeff(u_exact);
+   // Given advection vector:
+   VectorFunctionCoefficient advection(dim, advection_function);
+
+   // 6. Define the different forms and gridfunctions.
+
+   // We apply static condensation to the system
+   //
+   // [ A   B ] [  u   ] = [ F ]
+   // [ C   D ] [ ubar ]   [ H ]
+   //
+   // Eliminating u we find the global system
+   //
+   // S ubar = G
+   //
+   // where S = - C A^{-1} B + D and G = -C A^{-1} F + H.
+   // Having solved this system for ubar, we can compute u from
+   //
+   // u = A^{-1} (F - B ubar)
+
+   // Set up the linear form fform(.) which corresponds to the right-hand
+   // side of the linear system, which in this case is (f, phi_i) and
+   // phi_i are the basis functions in the finite element Uh_space.
+   FunctionCoefficient fcoeff(f_rhs);
+   LinearForm *fform(new LinearForm);
+   fform->AddDomainIntegrator(new DomainLFIntegrator(fcoeff));
+
+   // Set up the linear form gform(.) which corresponds to the right-hand
+   // side of the linear system, which in this case is <g, bar_phi_i>_{Gamma_N} and
+   // bar_phi_i are the basis functions in the finite element Uhbar_space.
+   LinearForm *gform(new LinearForm);
+   gform->AddSktBoundaryNeumannIntegrator(new HDGInflowLFIntegrator(ucoeff,
+                                                                    advection));
+
+   // Set up the bilinear form for the whole system. HDGBilinearForm2 can compute
+   // the Schur complement locally for a 2x2 problem.
+   HDGBilinearForm *AVarf(new HDGBilinearForm(Uh_space, Uhbar_space));
+   AVarf->AddHDGDomainIntegrator(
+      new HDGDomainIntegratorAdvection(mu, advection));
+   AVarf->AddHDGFaceIntegrator(
+      new HDGFaceIntegratorAdvection(advection));
+
+   GridFunction ubar(Uhbar_space);
+   GridFunction u(Uh_space);
+
+   for (int ref_levels = 0; ref_levels < total_ref_levels; ref_levels++)
+   {
+      // 7. Define the right hand side vectors
+      int dimUh   = Uh_space->GetVSize();
+      int dimUhbar = Uhbar_space->GetVSize();
+
+      std::cout << "***********************************************************\n";
+      std::cout << "dim(Uh) = " << dimUh << "\n";
+      std::cout << "dim(Uhbar) = " << dimUhbar << "\n";
+      std::cout << "***********************************************************\n";
+
+      Vector rhs_F(dimUh);
+      Vector rhs_G(dimUhbar);
+      Vector UBAR(dimUhbar);
+
+      // 8. Assemble the RHS and the bilinear forms
+      fform->Update(Uh_space, rhs_F, 0);
+      fform->Assemble(); // This is a vector
+
+      gform->Update(Uhbar_space, rhs_G, 0);
+      gform->Assemble(); // This is a vector
+
+      // Compute and Finalize the Schur complement
+      GridFunction *F = new GridFunction(Uh_space, rhs_F);
+      AVarf->AssembleSC(F, memA, memB);
+      AVarf->Finalize();
+
+      SparseMatrix *SC = AVarf->SpMatSC();
+
+      Vector *rhs_SC = AVarf->VectorSC();
+
+      // AVarf->VectorSC() provides -C*A^{-1} F, but the RHS for the
+      // Schur complement is  G - C*A^{-1} F
+      *rhs_SC += rhs_G;
+
+      // 9. Solve the Schur complement system
+      const int maxIter(1000);
+      const real_t rtol(1.e-15);
+      const real_t atol(0.0);
+      const int PrintLevel = -1;
+      GSSmoother M(*SC, 1, 1);
+      BiCGSTABSolver solver;
+      solver.SetAbsTol(atol);
+      solver.SetRelTol(rtol);
+      solver.SetMaxIter(maxIter);
+      solver.SetOperator(*SC);
+      solver.SetPrintLevel(PrintLevel);
+      solver.SetPreconditioner(M);
+      ubar = 0.0;
+      chrono.Clear();
+      chrono.Start();
+      solver.Mult(*rhs_SC, ubar);
+      chrono.Stop();
+
+      if (solver.GetConverged())
+         std::cout << "Iterative method converged in "
+                   << solver.GetNumIterations()
+                   << " iterations with a residual norm of "
+                   << solver.GetFinalNorm() << ".\n";
+      else
+         std::cout << "Iterative method did not converge in "
+                   << solver.GetNumIterations()
+                   << " iterations. Residual norm is "
+                   << solver.GetFinalNorm() << ".\n";
+
+      std::cout << "Iterative method solver took "
+                << chrono.RealTime() << "s. \n";
+      iterativeMethodIts[ref_levels] = solver.GetNumIterations();
+
+      // Delete the SC matrix to save memory
+      SC = NULL;
+
+      // 10. Reconstruction
+      // Reconstruct the solution u from the facet solution ubar
+      AVarf->Reconstruct(F, &ubar, &u);
+
+      // 11. Compute the discretization error
+      const int order_quad = max(2, 2*order+2);
+      const IntegrationRule *irs[Geometry::NumGeom];
+      for (int i=0; i < Geometry::NumGeom; ++i)
+      {
+         irs[i] = &(IntRules.Get(i, order_quad));
+      }
+
+      const real_t err_u  = u.ComputeL2Error(ucoeff, irs);
+      l2errors(ref_levels) = fabs(
+                                err_u); // fabs() to avoid negative values that ComputeL2Error can create
+
+      // 12. Save the mesh and the solution.
+      if (save)
+      {
+         ofstream mesh_ofs("refined.mesh");
+         mesh_ofs.precision(8);
+         mesh->Print(mesh_ofs);
+
+         ofstream u_ofs("sol_u.gf");
+         u_ofs.precision(8);
+         u.Save(u_ofs);
+
+         ofstream ubar_ofs("sol_lambda.gf");
+         ubar_ofs.precision(8);
+         ubar.Save(ubar_ofs);
+      }
+
+      // 13. Send the solution by socket to a GLVis server.
+      if (visualization)
+      {
+         char vishost[] = "localhost";
+         int  visport   = 19916;
+         socketstream u_sock(vishost, visport);
+         u_sock.precision(8);
+         u_sock << "solution\n" << *mesh << u << flush;
+      }
+
+      // 14. Refine the mesh to increase the resolution and update the spaces and the forms.
+      mesh->UniformRefinement();
+
+      Uh_space->Update(0);
+      Uhbar_space->Update(0);
+
+      AVarf->Update();
+
+      u.Update();
+      ubar.Update();
+
+      delete F;
+   }
+
+   // 15. Print the results and compute the rates
+   std::cout << "\n\n---------------------------------\n";
+   std::cout << "level  l2errors  order iterations\n";
+   std::cout << "---------------------------------\n";
+   for (int ref_levels = 0; ref_levels < total_ref_levels; ref_levels++)
+   {
+      if (ref_levels == 0)
+      {
+         std::cout << "  " << ref_levels << "   "
+                   << std::setprecision(2) << std::scientific
+                   << l2errors(ref_levels)
+                   << "  " << "-    " << "    "
+                   << iterativeMethodIts[ref_levels] << std::endl;
+      }
+      else
+      {
+         const real_t conv_order = log(l2errors(ref_levels)/l2errors(ref_levels-1))/log(
+                                      0.5);
+         std::cout << "  " << ref_levels << "   "
+                   << std::setprecision(2) << std::scientific
+                   << l2errors(ref_levels)
+                   << "  " << std::setprecision(4) << std::fixed
+                   << conv_order << "   "
+                   << iterativeMethodIts[ref_levels] << std::endl;
+      }
+   }
+   std::cout << "\n\n";
+
+   // 16. Free the used memory.
+   delete mesh;
+   delete Uh_fec;
+   delete Uhbar_fec;
+   delete Uh_space;
+   delete Uhbar_space;
+   delete fform;
+   delete gform;
+   delete AVarf;
+
+   std::cout << "Done." << std::endl ;
+
+   return 0;
+}
+//---------------------------------------------------------------------
+// Exact solution
+real_t u_exact(const Vector &x)
+{
+   real_t ue = 0.0;
+   const real_t xx = x(0);
+   const real_t yy = x(1);
+   if (dim == 2)
+   {
+      ue = 1.0 + sin(0.125 * M_PI * (1.0+xx) * (1.0+yy) * (1.0+yy));
+   }
+   else if (dim == 3)
+   {
+      const real_t zz = x(2);
+      ue = 1.0 + sin(0.125 * M_PI * (1.0+xx) * (1.0+yy) * (1.0+zz));
+   }
+
+   return ue;
+}
+//---------------------------------------------------------------------
+// Rhs function
+real_t f_rhs(const Vector &x)
+{
+   real_t rhs = 0.0;
+   const real_t ax = 0.8;
+   const real_t ay = 0.6;
+   const real_t mu = 1.0;
+   const real_t xx = x(0);
+   const real_t yy = x(1);
+
+   if (dim == 2)
+   {
+      const real_t uu = 1.0 + sin(0.125 * M_PI * (1.0+xx) * (1.0+yy) * (1.0+yy));
+      const real_t dudx = 0.125 * M_PI * (1.0+yy) * (1.0+yy)
+                          * cos(0.125 * M_PI * (1.0+xx) * (1.0+yy) * (1.0+yy));
+      const real_t dudy =  0.25 * M_PI * (1.0+xx) * (1.0+yy)
+                           * cos(0.125 * M_PI * (1.0+xx) * (1.0+yy) * (1.0+yy));
+
+      rhs = mu * uu + ax * dudx + ay * dudy;
+   }
+
+   if (dim == 3)
+   {
+      const real_t az = 0.7;
+      const real_t zz = x(2);
+      const real_t uu = 1.0 + sin(0.125 * M_PI * (1.0+xx) * (1.0+yy) * (1.0+zz));
+      const real_t dudx = 0.125 * M_PI * (1.0+yy) * (1.0+zz)
+                          * cos(0.125 * M_PI * (1.0+xx) * (1.0+yy) * (1.0+zz));
+
+      const real_t dudy = 0.125 * M_PI * (1.0+xx) * (1.0+zz)
+                          * cos(0.125 * M_PI * (1.0+xx) * (1.0+yy) * (1.0+zz));
+
+      const real_t dudz = 0.125 * M_PI * (1.0+xx) * (1.0+yy)
+                          * cos(0.125 * M_PI * (1.0+xx) * (1.0+yy) * (1.0+zz));
+
+      rhs = mu* uu + ax * dudx + ay * dudy + az * dudz;
+   }
+
+   return rhs;
+}
+//---------------------------------------------------------------------
+// Advection vector
+void advection_function(const Vector &x, Vector &v)
+{
+   if (dim == 3)
+   {
+      v(0) = 0.8;
+      v(1) = 0.6;
+      v(2) = 0.7;
+   }
+   else if (dim == 2)
+   {
+      v(0) = 0.8;
+      v(1) = 0.6;
+   }
+}

miniapps/hdg/hdg_advectionp.cpp

@@ -0,0 +1,606 @@
+//                        MFEM EDG/HDG example
+//
+// Compile with: make advectionp
+//
+// Sample runs:  mpirun -np 1 hdg_advectionp -o 1 -r 4 -tr 1 -no-vis // test scalability
+//            mpirun -np 2 hdg_advectionp -o 1 -r 4 -tr 1 -no-vis // test scalability
+//            mpirun -np 4 hdg_advectionp -o 1 -r 4 -tr 1 -no-vis // test scalability
+//            mpirun -np 4 hdg_advectionp -o 1 -r 0 -tr 3 -no-vis // test conv rates
+//            mpirun -np 4 hdg_advectionp -o 1 -r 2 -tr 1 -no-vis // test scalability
+//            mpirun -np 2 hdg_advectionp -o 5 -r 4 -tr 1 -no-vis // test conv rates
+//            mpirun -np 2 hdg_advectionp -o 5 -r 5 -tr 1 -m ../data/inline-tri.mesh -no-vis // test
+//
+// Description:  This example code demonstrates the use of MFEM to define a
+//               finite element discretization of the advection-reaction problem
+//               mu u + a.grad(u) = f with inhomogeneous Neumann boundary conditions.
+//               Specifically, we discretize using a HDG space of the
+//               specified order.
+//
+// The weak form is: seek (u,ubar) such that for all (v, vbar)
+//
+// \mu (u,v)   + (v, a.grad(u))  - < 1, [zeta a.n u v] >  + < ubar, [zeta a.n v] > = (f, w)
+// < ubar, [zeta a.n v] >  +
+//       < 1, [zeta a.n ubar vbar] > + < 1, [(1-zeta) a.n ubar vbar >_{\Gamma_N} = < g, vbar >
+//
+// where (.,.) is the d-dimensional L2 product, <.,.> is the d-1 dimensional L2 product.
+//
+// The discretization is based on the paper:
+//
+// G. N. Wells, Analysis of an interface stabilized finite element method: the advection-diffusion-reaction equation, SIAM J. Numer. Anal., 2011, 49:1, 87--109.
+//
+// Contributed by: T. Horvath, S. Rhebergen, A. Sivas
+//             University of Waterloo
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include <algorithm>
+#include "HDGBilinearForm.hpp"
+#include "hdg_integrators.hpp"
+
+using namespace std;
+using namespace mfem;
+
+//---------------------------------------------------------------------
+// Exact solution and r.h.s.. See below for implementation.
+real_t u_exact(const Vector &x);
+real_t f_rhs  (const Vector &x);
+void advection_function(const Vector &x, Vector &v);
+int dim;
+//---------------------------------------------------------------------
+int main(int argc, char *argv[])
+{
+   StopWatch chrono;
+
+   // 1. Initialize MPI.
+   // 1. Initialize MPI.
+   Mpi::Init(argc, argv);
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();
+
+   real_t assemblyTime, solveTime, reconstructTime;
+   real_t GassemblyTime, GsolveTime, GreconstructTime;
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../../data/inline-tri.mesh";
+   int order          = 1;
+   int initial_ref_levels = 0;
+   int total_ref_levels  = 2;
+   bool visualization   = true;
+   bool save          = true;
+   bool hdg = true;
+   real_t memA = 0.0;
+   real_t memB = 0.0;
+   bool petsc = false;
+   bool verbose = (myid == 0);
+   const char *petscrc_file = "";
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree > 1).");
+   args.AddOption(&initial_ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly for the initial calculation.");
+   args.AddOption(&total_ref_levels, "-tr", "--refine",
+                  "Number of times to refine the mesh uniformly.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&save, "-save", "--save-files", "-no-save",
+                  "--no-save-files",
+                  "Enable or disable file saving.");
+   args.AddOption(&hdg, "-hdg", "--hybrid", "-edg",
+                  "--embedded",
+                  "HDG / EDG option.");
+   args.AddOption(&memA, "-memA", "--memoryA",
+                  "Storage of A.");
+   args.AddOption(&memB, "-memB", "--memoryB",
+                  "Storage of B.");
+   args.AddOption(&petsc, "-petsc", "--use-petsc",
+                  "-no-petsc", "--no-use-petsc",
+                  "Enable or disable SC solver.");
+   args.AddOption(&petscrc_file, "-petscopts", "--petscopts",
+                  "PetscOptions file to use.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+#ifdef MFEM_USE_PETSC
+   // We initialize PETSc
+   MFEMInitializePetsc(NULL,NULL,petscrc_file,NULL);
+#endif
+
+#ifndef MFEM_USE_PETSC
+   if (petsc)
+   {
+      std::cout << "MFEM does not use PETSc. Change the solver to hypre" << std::endl
+                << std::flush;
+      petsc = false;
+   }
+#endif
+
+   if (order < 1)
+   {
+      cout << "Polynomial order should be > 0. Changing to order 1.";
+      order = 1;
+   }
+
+   if (memB > memA)
+   {
+      std::cout << "memB cannot be more than memA. Resetting to be equal" << std::endl
+                << std::flush;
+      memA = memB;
+   }
+   if (memA > 1.0)
+   {
+      std::cout << "memA cannot be more than 1. Resetting to 1" << std::endl <<
+                std::flush;
+      memA = 1.0;
+   }
+   else if (memA < 0.0)
+   {
+      std::cout << "memA cannot be less than 0. Resetting to 0." << std::endl <<
+                std::flush;
+      memA = 0.0;
+   }
+   if (memB > 1.0)
+   {
+      std::cout << "memB cannot be more than 1. Resetting to 1" << std::endl <<
+                std::flush;
+      memB = 1.0;
+   }
+   else if (memB < 0.0)
+   {
+      std::cout << "memB cannot be less than 0. Resetting to 0." << std::endl <<
+                std::flush;
+      memB = 0.0;
+   }
+
+   // 3. Read the mesh from the given mesh file. Refine it up to the initial_ref_levels.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   dim = mesh->Dimension();
+
+   if (mesh->Nonconforming())
+   {
+      if (verbose)
+      {
+         cout << "The current implementation does not support Nonconforming meshes. Terminating"
+              << endl << flush;
+      }
+#ifdef MFEM_USE_PETSC
+      MFEMFinalizePetsc();
+#endif
+      MPI_Finalize();
+      return 1;
+   }
+
+
+   for (int ii=0; ii<initial_ref_levels; ii++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // Define a parallel mesh. The serial mesh can be deleted.
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+
+   // 4. Define the vectors that will contain the errors and the iteration count at every refinement level
+   Vector l2errors(total_ref_levels);
+   Array<int> iterativeMethodIts(total_ref_levels);
+
+   // Define parallel finite element spaces on the parallel mesh.
+   // Uh_space is the DG space on elements
+   // ubar_space is the DG space on faces
+
+   if (order < 1)
+   {
+      cout << "Polynomial order should be > 0. Changing to order 1.";
+      order = 1;
+   }
+
+   // 5. Define the finite element spaces on the mesh.
+   FiniteElementCollection *Uh_fec(new DG_FECollection(order, dim));
+   FiniteElementCollection *Uhbar_fec = NULL;
+   if (hdg)
+   {
+      Uhbar_fec = new DG_Interface_FECollection(order, dim);
+   }
+   else
+   {
+      Uhbar_fec = new H1_Trace_FECollection(order, dim);
+   }
+
+   ParFiniteElementSpace *Uh_space   = new ParFiniteElementSpace(pmesh, Uh_fec);
+   ParFiniteElementSpace *Uhbar_space = new ParFiniteElementSpace(pmesh,
+                                                                  Uhbar_fec);
+
+   // 6. Define the coefficients
+   FunctionCoefficient fcoeff(f_rhs);
+   FunctionCoefficient ucoeff(u_exact);
+   ConstantCoefficient mu(1.0); // reaction constant
+
+   // Given advection vector:
+   VectorFunctionCoefficient advection(dim, advection_function);
+
+   // 7. Define the different forms and gridfunctions.
+
+   // Set up the linear form fform(.) which corresponds to the right-hand
+   // side of the linear system, which in this case is (f, phi_i) and
+   // phi_i are the basis functions in the finite element Uh_space.
+
+   ParLinearForm *fform = new ParLinearForm(Uh_space);
+   fform->AddDomainIntegrator(new DomainLFIntegrator(fcoeff));
+
+   // Set up the linear form gform(.) which corresponds to the right-hand
+   // side of the linear system, which in this case is <g, bar_phi_i>_{Gamma_N} and
+   // bar_phi_i are the basis functions in the finite element Uhbar_space.
+
+   ParLinearForm *gform = new ParLinearForm(Uhbar_space);
+   gform->AddSktBoundaryNeumannIntegrator(new HDGInflowLFIntegrator(ucoeff,
+                                                                    advection));
+
+   // We apply static condensation to the system
+   //
+   // [ A   B ] [  u   ] = [ F ]
+   // [ C   D ] [ ubar ]   [ H ]
+   //
+   // Eliminating u we find the global system
+   //
+   // S ubar = G
+   //
+   // where S = - C A^{-1} B + D and G = -C A^{-1} F + H.
+   // Having solved this system for ubar, we can compute u from
+   //
+   // u = A^{-1} (F - B ubar)
+
+   // Set up the bilinear form for the whole system. ParHDGBilinearForm2 can compute
+   // the Schur complement locally for a 2x2 problem.
+   HDGBilinearForm *AVarf(new HDGBilinearForm(Uh_space, Uhbar_space, true));
+   AVarf->AddHDGDomainIntegrator(
+      new HDGDomainIntegratorAdvection(mu, advection));
+   AVarf->AddHDGFaceIntegrator(
+      new HDGFaceIntegratorAdvection(advection));
+
+   ParGridFunction u(Uh_space);
+   ParGridFunction ubar(Uhbar_space);
+
+   for (int ref_levels = 0; ref_levels < total_ref_levels; ref_levels++)
+   {
+      // 8. Define the right hand side vectors
+      HYPRE_Int dimUh   = Uh_space->GlobalTrueVSize();
+      HYPRE_Int dimUhbar = Uhbar_space->GlobalTrueVSize();
+
+      if (verbose)
+      {
+         std::cout << "****************************************************\n";
+         std::cout << "dim(Uh)      = " << dimUh << "\n";
+         std::cout << "dim(Uhbar)   = " << dimUhbar << "\n";
+         std::cout << "dim(Uh+Uhbar) = " << dimUh + dimUhbar << "\n";
+         std::cout << "****************************************************\n";
+      }
+
+      ubar = 0.0;
+      HypreParVector *UBAR = ubar.ParallelProject();
+
+      // 10. Assemble the RHS and the bilinear forms
+      Vector rhs_F(Uh_space->GlobalVSize());
+      Vector rhs_H(Uhbar_space->GlobalVSize());
+
+      // Linear forms
+      fform->Update(Uh_space);
+      fform->Assemble();
+      gform->Update(Uhbar_space);
+      gform->Assemble();
+
+      HypreParVector *trueF;
+      trueF = fform->ParallelAssemble();
+
+      HypreParVector *trueG;
+      trueG = gform->ParallelAssemble();
+      // Create a ParGridFunction from the right hand side.
+      ParGridFunction *F = new ParGridFunction(Uh_space, trueF);
+
+      chrono.Clear();
+      chrono.Start();
+      AVarf->AssembleSC(F, memA, memB);
+      chrono.Stop();
+      AVarf->Finalize();
+
+      assemblyTime = chrono.RealTime();
+
+      HypreParMatrix *SC = AVarf->ParallelAssembleSC();
+
+      HypreParVector *rhs_SC = AVarf->ParallelVectorSC();
+      // AVarf->ParallelVectorSC() provides -C*A^{-1} F, but the RHS for the
+      // Schur complement is  G - C*A^{-1} F
+      rhs_SC->Add(1.0, *trueG);
+
+      real_t tol = 1.0e-12;
+      int maxIter = 1000;
+      int PrintLevel = -1;
+
+      // 12. Solve the linear system
+      if (petsc)
+      {
+#ifdef MFEM_USE_PETSC
+         // Solver using PETSc
+         //=======================
+         PetscLinearSolver *petsc_solver;
+         PetscPreconditioner *petsc_precon= NULL;
+         petsc_solver = new PetscLinearSolver(MPI_COMM_WORLD, "solver_");
+         petsc_precon = new PetscPreconditioner(MPI_COMM_WORLD,*SC,"solver_");
+         petsc_solver->SetOperator(*SC);
+         petsc_solver->SetPreconditioner(*petsc_precon);
+         petsc_solver->SetTol(tol);
+         petsc_solver->SetAbsTol(0.0);
+         petsc_solver->SetMaxIter(maxIter);
+         petsc_solver->SetPrintLevel(PrintLevel);
+         chrono.Clear();
+         chrono.Start();
+         petsc_solver->Mult(*rhs_SC, *UBAR);
+         chrono.Stop();
+
+         if (verbose)
+         {
+            if (petsc_solver->GetConverged())
+               std::cout << "Solver converged in " << petsc_solver->GetNumIterations()
+                         << " iterations with a residual norm of " << petsc_solver->GetFinalNorm() <<
+                         ".\n";
+            else
+               std::cout << "Solver did not converge in " << petsc_solver->GetNumIterations()
+                         << " iterations. Residual norm is " << petsc_solver->GetFinalNorm() << ".\n";
+
+            std::cout << "Solver solver took " << chrono.RealTime() << "s. \n";
+         }
+         delete petsc_solver;
+         delete petsc_precon;
+#endif
+      }
+      else
+      {
+         HypreSolver *pdiag = new HypreDiagScale(*SC);
+         HypreGMRES *itsolver = new HypreGMRES(*SC);
+         itsolver->SetTol(tol);
+         itsolver->SetMaxIter(maxIter);
+         itsolver->SetPrintLevel(PrintLevel);
+         itsolver->SetPreconditioner(*pdiag);
+         chrono.Clear();
+         chrono.Start();
+         itsolver->Mult(*rhs_SC, *UBAR);
+         chrono.Stop();
+
+         int numIterations = 0;
+         itsolver->GetNumIterations(numIterations);
+         if (verbose)
+         {
+            std::cout << "\nIterative method converged in "
+                      << numIterations << ".\n";
+
+            iterativeMethodIts[ref_levels] = numIterations;
+            std::cout << "Iterative solver took " << chrono.RealTime() << "s. \n";
+         }
+      }
+
+      // Delete the SC matrix to save memory
+      SC = NULL;
+      solveTime = chrono.RealTime();
+
+      ubar = ParGridFunction(Uhbar_space, UBAR);
+
+      chrono.Clear();
+      chrono.Start();
+      AVarf->Reconstruct(F, &ubar, &u);
+      chrono.Stop();
+
+      reconstructTime = chrono.RealTime();
+
+      const int order_quad = max(2, 2*order+2);
+      const IntegrationRule *irs[Geometry::NumGeom];
+      for (int i=0; i < Geometry::NumGeom; ++i)
+      {
+         irs[i] = &(IntRules.Get(i, order_quad));
+      }
+
+      const real_t err_u  = u.ComputeL2Error(ucoeff, irs);
+      if (verbose)
+      {
+         std::cout << "\nL2 error " << err_u << ".\n";
+      }
+      l2errors(ref_levels) = fabs(
+                                err_u); // fabs() to avoid negative values that ComputeL2Error can create
+
+      // 14. Save the mesh and the solution.
+      if (save)
+      {
+         ostringstream mesh_name, u_name, ubar_name;
+         mesh_name << "mesh." << setfill('0') << setw(6) << myid;
+         u_name << "sol_u." << setfill('0') << setw(6) << myid;
+         ubar_name << "sol_ubar." << setfill('0') << setw(6) << myid;
+
+         ofstream mesh_ofs(mesh_name.str().c_str());
+         mesh_ofs.precision(8);
+         pmesh->Print(mesh_ofs);
+
+         ofstream u_ofs(u_name.str().c_str());
+         u_ofs.precision(8);
+         u.Save(u_ofs);
+      }
+
+      // 15. Send the solution by socket to a GLVis server.
+      if (visualization)
+      {
+         char vishost[] = "localhost";
+         int  visport   = 19916;
+         socketstream u_sock(vishost, visport);
+         u_sock << "parallel " << num_procs << " " << myid << "\n";
+         u_sock.precision(8);
+         u_sock << "solution\n" << *pmesh << u << "window_title 'Velocity'"
+                << endl;
+         // Make sure all ranks have sent their 'u' solution before initiating
+         // another set of GLVis connections (one from each rank):
+         MPI_Barrier(pmesh->GetComm());
+      }
+
+      // 16. Refine the mesh to increase the resolution and update the spaces and the forms. Print the runtimes
+      pmesh->UniformRefinement();
+
+      Uh_space->Update(0);
+      Uhbar_space->Update(0);
+
+      AVarf->Update();
+
+      u.Update();
+      ubar.Update();
+      MPI_Reduce(&assemblyTime,&GassemblyTime,1,MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD);
+      MPI_Reduce(&solveTime,&GsolveTime,1,MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD);
+      MPI_Reduce(&reconstructTime,&GreconstructTime,1,MPI_DOUBLE,MPI_MAX,0,
+                 MPI_COMM_WORLD);
+
+      if (verbose)
+      {
+         printf("\t Assembly time   = %.2f\n",GassemblyTime);
+         printf("\t Solve time      = %.2f\n",GsolveTime);
+         printf("\t Reconstruct time = %.2f\n",GreconstructTime);
+      }
+
+      delete F;
+   }
+
+   // 17. Print the results
+   if (verbose)
+   {
+      std::cout << "\n\n---------------------------------\n";
+      std::cout << "level  l2errors  order iterations\n";
+      std::cout << "---------------------------------\n";
+      for (int ref_levels = 0; ref_levels < total_ref_levels; ref_levels++)
+      {
+         if (ref_levels == 0)
+         {
+            std::cout << "  " << ref_levels << "   "
+                      << std::setprecision(2) << std::scientific
+                      << l2errors(ref_levels)
+                      << "  " << "-    " << "    "
+                      << iterativeMethodIts[ref_levels] << std::endl;
+         }
+         else
+         {
+            const real_t conv_order = log(l2errors(ref_levels)/l2errors(ref_levels-1))
+                                      /log(0.5);
+            std::cout << "  " << ref_levels << "   "
+                      << std::setprecision(2) << std::scientific
+                      << l2errors(ref_levels)
+                      << "  " << std::setprecision(4) << std::fixed
+                      << conv_order << "   "
+                      << iterativeMethodIts[ref_levels] << std::endl;
+         }
+      }
+      std::cout << "\n\n";
+   }
+
+
+   // 19. Free the used memory.
+   delete pmesh;
+   delete Uh_fec;
+   delete Uhbar_fec;
+   delete Uh_space;
+   delete Uhbar_space;
+   delete fform;
+   delete gform;
+   delete AVarf;
+
+#ifdef MFEM_USE_PETSC
+   MFEMFinalizePetsc();
+#endif
+   MPI_Finalize();
+
+   return 0;
+}
+//---------------------------------------------------------------------
+// Exact solution
+real_t u_exact(const Vector &x)
+{
+   real_t ue = 0.0;
+   const real_t xx = x(0);
+   const real_t yy = x(1);
+   if (dim == 2)
+   {
+      ue = 1.0 + sin(0.125 * M_PI * (1.0+xx) * (1.0+yy) * (1.0+yy));
+   }
+   else if (dim == 3)
+   {
+      const real_t zz = x(2);
+      ue = 1.0 + sin(0.125 * M_PI * (1.0+xx) * (1.0+yy) * (1.0+zz));
+   }
+
+   return ue;
+}
+//---------------------------------------------------------------------
+// Rhs function
+real_t f_rhs(const Vector &x)
+{
+   real_t rhs = 0.0;
+   const real_t ax = 0.8;
+   const real_t ay = 0.6;
+   const real_t mu = 1.0;
+   const real_t xx = x(0);
+   const real_t yy = x(1);
+
+   if (dim == 2)
+   {
+      const real_t uu = 1.0 + sin(0.125 * M_PI * (1.0+xx) * (1.0+yy) * (1.0+yy));
+      const real_t dudx = 0.125 * M_PI * (1.0+yy) * (1.0+yy)
+                          * cos(0.125 * M_PI * (1.0+xx) * (1.0+yy) * (1.0+yy));
+      const real_t dudy =  0.25 * M_PI * (1.0+xx) * (1.0+yy)
+                           * cos(0.125 * M_PI * (1.0+xx) * (1.0+yy) * (1.0+yy));
+
+      rhs = mu * uu + ax * dudx + ay * dudy;
+   }
+
+   if (dim == 3)
+   {
+      const real_t az = 0.7;
+      const real_t zz = x(2);
+      const real_t uu = 1.0 + sin(0.125 * M_PI * (1.0+xx) * (1.0+yy) * (1.0+zz));
+      const real_t dudx = 0.125 * M_PI * (1.0+yy) * (1.0+zz)
+                          * cos(0.125 * M_PI * (1.0+xx) * (1.0+yy) * (1.0+zz));
+
+      const real_t dudy = 0.125 * M_PI * (1.0+xx) * (1.0+zz)
+                          * cos(0.125 * M_PI * (1.0+xx) * (1.0+yy) * (1.0+zz));
+
+      const real_t dudz = 0.125 * M_PI * (1.0+xx) * (1.0+yy)
+                          * cos(0.125 * M_PI * (1.0+xx) * (1.0+yy) * (1.0+zz));
+
+      rhs = mu * uu + ax * dudx + ay * dudy + az * dudz;
+   }
+
+   return rhs;
+}
+//---------------------------------------------------------------------
+// Advection vector
+void advection_function(const Vector &x, Vector &v)
+{
+   if (dim == 3)
+   {
+      v(0) = 0.8;
+      v(1) = 0.6;
+      v(2) = 0.7;
+   }
+   else if (dim == 2)
+   {
+      v(0) = 0.8;
+      v(1) = 0.6;
+   }
+}

miniapps/hdg/hdg_integrators.cpp

@@ -0,0 +1,615 @@
+// 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.
+//
+// Implementation of Bilinear Form Integrators
+//
+// Contributed by: T. Horvath: Oakland University
+//                 S. Rhebergen, A. Sivas: University of Waterloo
+
+#include "../../fem/fem.hpp"
+#include <cmath>
+#include <algorithm>
+#include "hdg_integrators.hpp"
+
+using namespace std;
+
+namespace mfem
+{
+void HDGDomainIntegratorAdvection::AssembleElementMatrix(
+   const FiniteElement &fe_u,
+   ElementTransformation &Trans,
+   DenseMatrix &elmat)
+{
+   int ndof_u = fe_u.GetDof();
+   int dim  = fe_u.GetDim();
+   int spaceDim = Trans.GetSpaceDim();
+   bool square = (dim == spaceDim);
+
+   Vector vec1; // for the convection integral
+   vec2.SetSize(dim);
+   BdFidxT.SetSize(ndof_u);
+
+   dshape.SetSize (ndof_u, dim); // for nabla \tilde{u}
+   gshape.SetSize (ndof_u, dim); // for nabla u
+   Jadj.SetSize (dim);  // for the Jacobian
+   shapeu.SetSize (ndof_u);    // shape of u
+
+   // setting the sizes of the local element matrices
+   elmat.SetSize(ndof_u, ndof_u);
+
+   // setting the order of integration
+   const IntegrationRule *ir = IntRule;
+   if (ir == NULL)
+   {
+      int order = 2 * fe_u.GetOrder() + 1;
+      ir = &IntRules.Get(fe_u.GetGeomType(), order);
+   }
+
+   elmat = 0.0;
+
+   // evaluate the advection vector at all integration point
+   avec->Eval(Adv_ir, Trans, *ir);
+
+   for (int i = 0; i < ir -> GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+
+      // shape functions
+      fe_u.CalcDShape (ip, dshape);
+      fe_u.CalcShape (ip, shapeu);
+
+      // calculate the Adjugate of the Jacobian
+      Trans.SetIntPoint (&ip);
+      CalcAdjugate(Trans.Jacobian(), Jadj);
+
+      real_t w = Trans.Weight();
+      w = ip.weight / (square ? w : w*w*w);
+      // AdjugateJacobian = / adj(J),         if J is square
+      //                    \ adj(J^t.J).J^t, otherwise
+
+      // Calculate the gradient of the function of the physical element
+      Mult (dshape, Jadj, gshape);
+
+      // get the advection at the current integration point
+      Adv_ir.GetColumnReference(i, vec1);
+      vec1 *= ip.weight; // so it will be (cu, nabla v)
+
+      // compute -(cu, nabla v)
+      Jadj.Mult(vec1, vec2);
+      dshape.Mult(vec2, BdFidxT);
+      AddMultVWt(shapeu, BdFidxT, elmat);
+
+      real_t massw = Trans.Weight() * ip.weight;
+
+      if (mass_coeff)
+      {
+         massw *= mass_coeff->Eval(Trans, ip);
+      }
+      AddMult_a_VVt(massw, shapeu, elmat);
+
+   }
+}
+//---------------------------------------------------------------------
+void HDGFaceIntegratorAdvection::AssembleFaceMatrixOneElement1and1FES(
+   const FiniteElement &fe_u,
+   const FiniteElement &face_fe,
+   FaceElementTransformations &Trans,
+   const int elem1or2,
+   const bool reconstruct_only,
+   DenseMatrix &elmat1,
+   DenseMatrix &elmat2,
+   DenseMatrix &elmat3,
+   DenseMatrix &elmat4)
+{
+   int dim, ndof, ndof_face;
+   real_t w;
+
+   dim = fe_u.GetDim();
+   ndof_face = face_fe.GetDof();
+
+   shape_face.SetSize(ndof_face);
+
+   normal.SetSize(dim);
+   normalJ.SetSize(dim);
+   invJ.SetSize(dim);
+   adv.SetSize(dim);
+
+   ndof = fe_u.GetDof();
+
+   shape.SetSize(ndof);
+   dshape.SetSize(ndof, dim);
+   dshape_normal.SetSize(ndof);
+
+   elmat1.SetSize(ndof);
+   elmat2.SetSize(ndof, ndof_face);
+   elmat3.SetSize(ndof_face, ndof);
+   elmat4.SetSize(ndof_face);
+
+   elmat1 = 0.0;
+   elmat2 = 0.0;
+   elmat3 = 0.0;
+   elmat4 = 0.0;
+
+   // Since we are using GetSharedFaceTransformations over the shard faces
+   // we can find the boundary elements by checking Trans.Elem2No
+   bool is_bdr = (Trans.Elem2No < 0);
+
+   const IntegrationRule *ir = IntRule;
+   if (ir == NULL)
+   {
+      // a simple choice for the integration order
+      int order;
+      order = 2*fe_u.GetOrder();
+
+      ir = &IntRules.Get(Trans.FaceGeom, order);
+   }
+
+   for (int p = 0; p < ir->GetNPoints(); p++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(p);
+      IntegrationPoint eip; // element integration point
+
+      Trans.Face->SetIntPoint(&ip);
+      face_fe.CalcShape(ip, shape_face);
+
+      if (dim == 1)
+      {
+         normal(0) = 2*eip.x - 1.0;
+      }
+      else
+      {
+         CalcOrtho(Trans.Face->Jacobian(), normal);
+      }
+
+      Trans.Loc1.Transform(ip, eip);
+      Trans.Elem1->SetIntPoint(&eip);
+
+      avec->Eval(adv, *Trans.Elem1, eip);
+      real_t an = adv * normal;
+      real_t an_L = an;
+
+      real_t zeta_R = 0.0, zeta_L = 0.0, zeta = 0.0;
+      if (an < 0.0)
+      {
+         zeta_L = 1.0;
+      }
+
+      if (elem1or2 == 1)
+      {
+         zeta = zeta_L;
+      }
+      else
+      {
+         Trans.Loc2.Transform(ip, eip);
+         Trans.Elem2->SetIntPoint(&eip);
+
+         avec->Eval(adv, *Trans.Elem2, eip);
+         an = adv * normal;
+         an *= -1.;
+
+         zeta_R = 1.0 - zeta_L;
+         zeta = zeta_R;
+      }
+
+      fe_u.CalcShape(eip, shape);
+
+      w = ip.weight;
+
+      for (int i = 0; i < ndof; i++)
+      {
+         for (int j = 0; j < ndof; j++)
+         {
+            // - < 1, [zeta a.n u v] >
+            elmat1(i, j) -= w * zeta * an * shape(i) * shape(j);
+         }
+
+         for (int j = 0; j < ndof_face; j++)
+         {
+            if (!reconstruct_only)
+            {
+               // - < ubar, [(1-zeta) a.n v] >
+               elmat3(j, i) -= w * an * (1.-zeta) * shape(i) * shape_face(j);
+            }
+
+            // + < ubar, [zeta a.n v] >
+            elmat2(i, j) += w * zeta * an * shape(i) * shape_face(j);
+
+         }
+      }
+      if (!reconstruct_only)
+      {
+
+         for (int i = 0; i < ndof_face; i++)
+            for (int j = 0; j < ndof_face; j++)
+            {
+               // - < 1, [zeta a.n ubar vbar] > + < 1, [(1-zeta) a.n ubar vbar >_{\Gamma_N}
+               if (!is_bdr)
+               {
+                  if (elem1or2 == 1)
+                  {
+                     elmat4(i, j) += -w * zeta_L * an_L * shape_face(i) * shape_face(j);
+                  }
+                  else
+                  {
+                     elmat4(i, j) += - w * (1.0 - zeta_L) * (-an_L) * shape_face(i) * shape_face(j);
+                  }
+               }
+               else
+               {
+                  elmat4(i, j) += -w * zeta_L * an * shape_face(i) * shape_face(j)
+                                  + w * (1.0 - zeta_L) * an * shape_face(i) * shape_face(j);
+               }
+            }
+      }
+
+   }
+}
+
+//---------------------------------------------------------------------
+void HDGInflowLFIntegrator::AssembleRHSElementVect(
+   const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
+{
+   mfem_error("Not implemented \n");
+}
+
+void HDGInflowLFIntegrator::AssembleRHSElementVect(
+   const FiniteElement &face_S, FaceElementTransformations &Trans,
+   Vector &favect)
+{
+   int dim, ndof_face;
+   real_t w, uin;
+
+   dim = face_S.GetDim(); // This is face dimension which is 1 less than
+   dim += 1;              // space dimension so add 1 to face dim to
+   // get the space dim.
+   n_L.SetSize(dim);
+   Vector adv(dim);
+
+   ndof_face = face_S.GetDof();
+
+   shape_f.SetSize(ndof_face);
+   favect.SetSize(ndof_face);
+   favect = 0.0;
+
+   if (Trans.Elem2No >= 0)
+   {
+      // Interior face, do nothing
+   }
+   else
+   {
+      // Boundary face
+      const IntegrationRule *ir = IntRule;
+      if (ir == NULL)
+      {
+         int order = 2 * face_S.GetOrder();
+         if (face_S.GetMapType() == FiniteElement::VALUE)
+         {
+            order += Trans.Face->OrderW();
+         }
+
+         ir = &IntRules.Get(Trans.FaceGeom, order);
+      }
+
+      for (int p = 0; p < ir->GetNPoints(); p++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(p);
+         face_S.CalcShape(ip, shape_f);
+
+         IntegrationPoint eip_L;
+         Trans.Loc1.Transform(ip, eip_L);
+         Trans.Face->SetIntPoint(&ip);
+
+         avec->Eval(adv, *Trans.Elem1, eip_L);
+         uin = u_in->Eval(*Trans.Elem1, eip_L);
+
+         if (dim == 1)
+         {
+            n_L(0) = 2*eip_L.x - 1.0;
+         }
+         else
+         {
+            CalcOrtho(Trans.Face->Jacobian(), n_L);
+         }
+
+         real_t an_L = adv * n_L;
+
+         real_t zeta_L = 0.0;
+         if (an_L < 0.0)
+         {
+            zeta_L = 1.0;
+         }
+
+         w = ip.weight;
+
+         real_t gg = -uin * an_L * zeta_L;
+         for (int i = 0; i < ndof_face; i++)
+         {
+            favect(i) += w * gg * shape_f(i);
+         }
+      }
+   }
+}
+//---------------------------------------------------------------------
+////////////////////////////////////////////////////////////////////////////////////////////////////
+void HDGDomainIntegratorDiffusion::AssembleElementMatrix2FES(
+   const FiniteElement &fe_q,
+   const FiniteElement &fe_u,
+   ElementTransformation &Trans,
+   DenseMatrix &elmat)
+{
+   // get the number of degrees of freedoms and the dimension of the problem
+   int ndof_u = fe_u.GetDof();
+   int ndof_q = fe_q.GetDof();
+   int dim  = fe_q.GetDim();
+   real_t norm;
+
+   int vdim = dim ;
+
+   // set the vector and matrix sizes
+   dshape.SetSize (ndof_u, dim); // for nabla u_reference
+   gshape.SetSize (ndof_u, dim); // for nabla u
+   Jadj.SetSize (dim);  // the Jacobian
+   divshape.SetSize (vdim*ndof_u);  // divergence of q
+   shape.SetSize (ndof_q);  // shape of q (and u)
+
+   // for vector diffusion the matrix is built up from partial matrices
+   partelmat.SetSize(ndof_q);
+
+   DenseMatrix local_A11, local_A12, local_A21;
+
+   // setting the sizes of the local element matrices
+   local_A11.SetSize(dim*ndof_q, dim*ndof_q);
+   local_A12.SetSize(vdim*ndof_q, ndof_u);
+   local_A21.SetSize(ndof_u, vdim*ndof_q);
+
+   elmat.SetSize(dim*ndof_q + ndof_u);
+
+   local_A11 = 0.0;
+   local_A12 = 0.0;
+   local_A21 = 0.0;
+   elmat  = 0.0;
+
+   // setting the order of integration
+   const IntegrationRule *ir = IntRule;
+   if (ir == NULL)
+   {
+      int order1 = 2 * fe_q.GetOrder();
+      int order2 = 2 * fe_q.GetOrder() + Trans.OrderW();
+      int order = max(order1, order2);
+      ir = &IntRules.Get(fe_u.GetGeomType(), order);
+   }
+
+   for (int i = 0; i < ir -> GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+
+      // compute the shape and the gradient values on the reference element
+      fe_u.CalcDShape (ip, dshape);
+      fe_q.CalcShape (ip, shape);
+
+      // calculate the adjugate of the Jacobian
+      Trans.SetIntPoint (&ip);
+      CalcAdjugate(Trans.Jacobian(), Jadj);
+
+      // Calculate the gradient of the function of the physical element
+      Mult (dshape, Jadj, gshape);
+
+      // the weight is the product of the integral weight and the
+      // determinant of the Jacobian
+      norm = ip.weight * Trans.Weight();
+      MultVVt(shape, partelmat);
+
+      real_t c = ip.weight;
+
+      // transform the the matrix to divergence vector
+      gshape.GradToDiv (divshape);
+
+      // mulitply by 1.0/nu
+      partelmat *= 1.0/nu->Eval(Trans, ip);
+
+      shape *= c;
+      // compute the (u, \div v) term
+      AddMultVWt (shape, divshape, local_A21);
+
+      // assemble -(q, v) from the partial matrices
+      partelmat *= norm*(-1.0);
+      for (int k = 0; k < vdim; k++)
+      {
+         local_A11.AddMatrix(partelmat, ndof_q*k, ndof_q*k);
+      }
+
+   }
+
+   local_A12.Transpose(local_A21);
+
+   int block_size1 = dim*ndof_q;
+
+   elmat.CopyMN(local_A11, 0, 0);
+   elmat.CopyMN(local_A12, 0, block_size1);
+   elmat.CopyMN(local_A21, block_size1, 0);
+}
+
+
+void HDGFaceIntegratorDiffusion::AssembleFaceMatrixOneElement2and1FES(
+   const FiniteElement &fe_q,
+   const FiniteElement &fe_u,
+   const FiniteElement &face_fe,
+   FaceElementTransformations &Trans,
+   const int elem1or2,
+   const bool reconstruct_only,
+   DenseMatrix &elmat1,
+   DenseMatrix &elmat2,
+   DenseMatrix &elmat3,
+   DenseMatrix &elmat4)
+{
+   // Get DoF from faces and the dimension
+   int ndof_face = face_fe.GetDof();
+   int ndof_q, ndof_u;
+   int dim = fe_q.GetDim();
+   int vdim = dim;
+   int order;
+
+   DenseMatrix shape1_n_mtx;
+
+   // set the dofs for u and q
+   ndof_u = fe_u.GetDof();
+   ndof_q = fe_q.GetDof();
+
+   DenseMatrix local_B1, local_A22, local_B2, local_C1, local_C2, local_D;
+
+   // set the shape functions, the normal and the advection
+   shapeu.SetSize(ndof_u);
+   shapeq.SetSize(ndof_q);
+   shape_face.SetSize(ndof_face);
+   normal.SetSize(dim);
+
+   // set the proper size for the matrices
+   local_B1.SetSize(vdim*ndof_q, ndof_face);
+   local_B1 = 0.0;
+   local_A22.SetSize(ndof_u, ndof_u);
+   local_A22 = 0.0;
+   local_B2.SetSize(ndof_u, ndof_face);
+   local_B2 = 0.0;
+   local_C1.SetSize(vdim*ndof_q, ndof_face);
+   local_C1 = 0.0;
+   local_C2.SetSize(ndof_u, ndof_face);
+   local_C2 = 0.0;
+   local_D.SetSize(ndof_face, ndof_face);
+   local_D = 0.0;
+
+   int sub_block_size1 = vdim*ndof_q;
+   int sub_block_size2 = ndof_u;
+
+   int block_size1 = sub_block_size1 + sub_block_size2;
+   int block_size2 = ndof_face;
+
+   elmat1.SetSize(block_size1);
+   elmat1 = 0.0;
+   elmat2.SetSize(block_size1, block_size2);
+   elmat2 = 0.0;
+   elmat3.SetSize(block_size2, block_size1);
+   elmat3 = 0.0;
+   elmat4.SetSize(block_size2);
+   elmat4 = 0.0;
+
+
+   shape1_n_mtx.SetSize(ndof_q,dim);
+   shape_dot_n.SetSize(ndof_q,dim);
+
+   // set the order of integration
+   // using the fact that q and u has the same order!
+   const IntegrationRule *ir = IntRule;
+   if (ir == NULL)
+   {
+      order = 2*max(max(fe_q.GetOrder(), fe_u.GetOrder()), face_fe.GetOrder());
+      order += 2;
+
+      // IntegrationRule depends on the Geometry of the face (pont, line, triangle, rectangular)
+      ir = &IntRules.Get(Trans.FaceGeom, order);
+   }
+
+   for (int p = 0; p < ir->GetNPoints(); p++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(p);
+      IntegrationPoint eip; // integration point on the element
+
+      // Trace finite element shape function
+      Trans.Face->SetIntPoint(&ip);
+      face_fe.CalcShape(ip, shape_face);
+
+      // calculate the normal at the integration point
+      if (dim == 1)
+      {
+         normal(0) = 2*eip.x - 1.0;
+      }
+      else
+      {
+         CalcOrtho(Trans.Face->Jacobian(), normal);
+      }
+
+      if (elem1or2 == 1)
+      {
+         // Side 1 finite element shape function
+         Trans.Loc1.Transform(ip, eip);
+      }
+      else
+      {
+         // Side 2 finite element shape function
+         Trans.Loc2.Transform(ip, eip);
+      }
+
+      fe_u.CalcShape(eip, shapeu);
+      fe_q.CalcShape(eip, shapeq);
+      MultVWt(shapeq, normal, shape_dot_n) ;
+
+
+      // set the coefficients for the different terms
+      // if the normal is involved, Trans.Face->Weight() is not required
+      real_t w1 = ip.weight*(-1.0);
+
+      if (elem1or2 == 2)
+      {
+         w1 *=-1.0;
+      }
+
+      real_t w2 = tauD*Trans.Face->Weight()* ip.weight;
+
+      real_t w3 = -w2;
+
+      // local_B1 = < \lambda,\nu v\cdot n>
+      for (int i = 0; i < vdim; i++)
+         for (int k = 0; k < ndof_q; k++)
+            for (int j = 0; j < ndof_face; j++)
+            {
+               local_B1(i*ndof_q + k, j) += shape_face(j) * shape_dot_n(k,i) * w1;
+            }
+
+      // local_A22 =  < \tau u, w>
+      // local_B2= -< tau \lambda, w>
+      // local_C2 = -< tau \lambda, w>
+      for (int i = 0; i < ndof_u; i++)
+      {
+         for (int j = 0; j < ndof_u; j++)
+         {
+            local_A22(i, j) += w2 * shapeu(i) * shapeu(j);
+         }
+
+         for (int j = 0; j < ndof_face; j++)
+         {
+            local_B2(i, j) += w3 * shapeu(i) * shape_face(j);
+         }
+      }
+
+      if (!reconstruct_only)
+      {
+         // local_D = < \tau \lambda, \mu>
+
+         AddMult_a_VVt(w2, shape_face, local_D);
+      }
+   }
+
+   local_C1.Transpose(local_B1);
+   local_C2.Transpose(local_B2);
+
+   elmat1.CopyMN(local_A22, sub_block_size1, sub_block_size1);
+
+   elmat2.CopyMN(local_B1, 0, 0);
+   elmat2.CopyMN(local_B2, sub_block_size1, 0);
+
+   elmat3.CopyMN(local_C1, 0, 0);
+   elmat3.CopyMN(local_C2, 0, sub_block_size1);
+
+   elmat4 = local_D;
+}
+
+
+}
+;

miniapps/hdg/hdg_integrators.hpp

@@ -0,0 +1,223 @@
+// 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.
+//
+// Implementation of Bilinear Form Integrators
+//
+// Contributed by: T. Horvath: Oakland University
+//                 S. Rhebergen, A. Sivas: University of Waterloo
+
+#ifndef MFEM_HDGINTEG
+#define MFEM_HDGINTEG
+
+#include "../../config/config.hpp"
+
+namespace mfem
+{
+
+//---------------------------------------------------------------------
+// Advection integrator: to compute all the domain based integrals
+//
+// The output is
+//
+// elemmat = mass_coeff (u,v) + (v, avec.grad(u))
+//
+// mass_coeff is the reaction coefficient
+// avec is the advection coefficient
+class HDGDomainIntegratorAdvection : public BilinearFormIntegrator
+{
+private:
+   Coefficient *mass_coeff;
+   VectorCoefficient *avec;
+
+   Vector shape1, shape2;
+   DenseMatrix shape1_n, shape2_n, partelmat;
+
+   Vector shapeq;
+   Vector shapeu;
+   Vector divshape, divshape_no_diffusion, vec2, BdFidxT;
+   DenseMatrix dshape;
+   DenseMatrix gshape;
+   DenseMatrix Jadj;
+   DenseMatrix Adv_ir;
+
+public:
+   HDGDomainIntegratorAdvection(Coefficient &mass, VectorCoefficient &_avec)
+      : mass_coeff(&mass), avec(&_avec)  {  }
+
+   using BilinearFormIntegrator::AssembleElementMatrix;
+   virtual void AssembleElementMatrix(const FiniteElement &fe_u,
+                                      ElementTransformation &Trans,
+                                      DenseMatrix &elmat1);
+};
+
+// Advection integrator to compute all the face based integrals
+//
+// The output is
+//
+// elemmat1 = - < 1, [zeta avec.n u v] >
+// elemmat2 = < ubar, [zeta avec.n v] >
+// elemmat3 = < ubar, [zeta avec.n v] >
+// elemmat4 = < 1, [zeta avec.n ubar vbar] > + < 1, [(1-zeta) avec.n ubar vbar >_{\Gamma_N}
+//
+// avec is the advection coefficient
+class HDGFaceIntegratorAdvection : public BilinearFormIntegrator
+{
+private:
+   VectorCoefficient *avec;
+
+   Vector shape, normal, normalJ, dshape_normal, shape_face, adv;
+   DenseMatrix shape1_n, invJ, dshape, shape1_n_nodiff;
+
+public:
+   HDGFaceIntegratorAdvection(VectorCoefficient &_avec)
+      : avec(&_avec)  {  }
+
+   using BilinearFormIntegrator::AssembleFaceMatrixOneElement1and1FES;
+   virtual void AssembleFaceMatrixOneElement1and1FES(const FiniteElement &fe_u,
+                                                     const FiniteElement &face_fe,
+                                                     FaceElementTransformations &Trans,
+                                                     const int elem1or2,
+                                                     const bool reconstruct_only,
+                                                     DenseMatrix &elmat1,
+                                                     DenseMatrix &elmat2,
+                                                     DenseMatrix &elmat3,
+                                                     DenseMatrix &elmat4);
+
+};
+
+//---------------------------------------------------------------------
+/** Boundary linear integrator for imposing inflow boundary
+    conditions. Given the inflow data u_in, the linear form assembles the
+    following integral on the boundary:
+
+    + < g, vbar >
+
+    where g = - u_in * a.n * zeta and vbar is the test function. */
+class HDGInflowLFIntegrator : public LinearFormIntegrator
+{
+protected:
+   Coefficient *u_in;
+   VectorCoefficient *avec;
+
+   // these are not thread-safe!
+   Vector shape_f, n_L;
+
+public:
+   HDGInflowLFIntegrator(Coefficient &_u, VectorCoefficient &_avec)
+   {
+      u_in = &_u;
+      avec = &_avec;
+   }
+
+   using LinearFormIntegrator::AssembleRHSElementVect;
+   virtual void AssembleRHSElementVect(const FiniteElement &el,
+                                       ElementTransformation &Tr,
+                                       Vector &elvect);
+
+   virtual void AssembleRHSElementVect(const FiniteElement &el,
+                                       FaceElementTransformations &Tr,
+                                       Vector &elvect);
+};
+
+//---------------------------------------------------------------------
+//---------------------------------------------------------------------
+
+// Diffusion integrator: to compute all the domain based integrals
+//
+// The output is
+//
+//         [local_A11 local_A12]
+// elmat = [local_A21    0.0   ]
+//
+// local_A11 = -(\nu^{-1} q, v)
+// local_A12 = (u, div(v))
+// local_A21 = (div(q), w)
+//
+// local_A21 = local_A12^T
+//
+// \nu is the constant diffusion coefficient
+class HDGDomainIntegratorDiffusion : public BilinearFormIntegrator
+{
+private:
+   ConstantCoefficient *nu;
+
+   Vector shape, divshape;
+   DenseMatrix partelmat, dshape, gshape, Jadj;
+
+public:
+   HDGDomainIntegratorDiffusion(ConstantCoefficient &_nu)
+      : nu(&_nu) { }
+
+   using BilinearFormIntegrator::AssembleElementMatrix2FES;
+   virtual void AssembleElementMatrix2FES(const FiniteElement &fe_q,
+                                          const FiniteElement &fe_u,
+                                          ElementTransformation &Trans,
+                                          DenseMatrix &elmat);
+
+};
+
+// Diffusion integrator to compute all the face based integrals
+//
+// The output is
+//
+//          [ 0.0    0.0   ]
+// elmat1 = [ 0.0 local_A22 ]  - the face based integral for matrix A
+//
+//          [ local_B1 ]
+// elmat2 = [ local_B2 ]  - the face based integral for matrix B
+//
+// elmat3 = [ local_C1  local_C2 ]  - the face based integral for matrix C
+//
+// elmat4 = local_D  - the face based integral for matrix D
+//
+// where
+// local_B1 = < \lambda,v\cdot n>
+// local_A22 = < \tau u, w>
+// local_B2 = -< tau \lambda, w>
+// local_C1 = < \lambda, v\cdot n>
+// local_C2 = -< \tau \lambda, w>
+// local_D = < \tau \lambda, \mu>
+//
+// q_diff_coeff is the constant diffusion coefficient
+// local_C1 = local_B1^T
+// local_C2 = local_B2^T
+class HDGFaceIntegratorDiffusion : public BilinearFormIntegrator
+{
+private:
+   real_t tauD;
+
+   Vector shapeu, shapeq, normal, shape_face;
+   DenseMatrix shape_dot_n;
+
+public:
+   HDGFaceIntegratorDiffusion(real_t a)
+   {
+      tauD = a;
+   }
+
+   using BilinearFormIntegrator::AssembleFaceMatrixOneElement2and1FES;
+   virtual void AssembleFaceMatrixOneElement2and1FES(const FiniteElement &fe_q,
+                                                     const FiniteElement &fe_u,
+                                                     const FiniteElement &face_fe,
+                                                     FaceElementTransformations &Trans,
+                                                     const int elem1or2,
+                                                     const bool reconstruct_only,
+                                                     DenseMatrix &elmat1,
+                                                     DenseMatrix &elmat2,
+                                                     DenseMatrix &elmat3,
+                                                     DenseMatrix &elmat4);
+
+};
+
+
+}
+
+#endif

miniapps/hdg/hdg_poisson.cpp

@@ -0,0 +1,675 @@
+//                        MFEM Example Hybridizable DG
+//
+// Compile with: make hdg_poisson
+//
+// Sample runs:  hdg_poisson -o 1 -r 1 -tr 4 -no-vis
+//            hdg_poisson -o 5 -r 1 -tr 4 -no-vis
+//            hdg_poisson -o 1 -r 4 -tr 1
+//            hdg_poisson -o 5 -r 4 -tr 1
+//            hdg_poisson -o 1 -r 1 -tr 4 -no-vis -m ../data/inline-tri.mesh
+//            hdg_poisson -o 5 -r 1 -tr 4 -no-vis -m ../data/inline-tri.mesh
+//            hdg_poisson -o 1 -r 5 -tr 1 -m ../data/inline-tri.mesh
+//            hdg_poisson -o 5 -r 5 -tr 1 -m ../data/inline-tri.mesh
+//
+// Description:  This example code solves the 2D/3D diffusion problem
+//                     -\nu Delta u = f
+//               with Dirichlet boundary conditions, using HDG discretization.
+//
+// The methods approximates the solution u, the diffusive flux q = -\nu \nabla u,
+// and the restriction of u to the faces, denoted by lambda.
+//
+// The weak form is: seek (q,u,\lambda) such that for all (v, w, \mu)
+//
+// -\nu^{-1}(q, v)       + (u, div(v))       - <\lambda, v \cdot n>      = 0
+//  (div(q), w)          + <\tau u, w>       - <\tau \lambda, w>         = (f, w)
+// -<[[q \cdot n]], \mu> - <[[\tau u]], \mu> + <[[(\tau \lambda]], \mu>  = 0
+//
+// where [[.]] is the jump operator, (.,.) is the d-dimensional L2 product,
+// <.,.> is the d-1 dimensional L2 product.
+//
+// The discretization is based on the paper:
+//
+// N.C. Nguyen, J. Peraire, B. Cockburn, An implicit high-order hybridizable
+// discontinuous Galerkin method for linear convection–diffusion equations,
+// J. Comput. Phys., 2009, 228:9, 3232--3254.
+//
+// Contributed by: T. Horvath, Oakland University
+//                 S. Rhebergen, A. Sivas, University of Waterloo
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include <algorithm>
+#include "HDGBilinearForm.hpp"
+#include "hdg_integrators.hpp"
+
+
+using namespace std;
+using namespace mfem;
+
+// Define the analytical solution and forcing terms / boundary conditions
+real_t uFun_ex(const Vector & x);
+void qFun_ex(const Vector & x, Vector & q);
+real_t fFun(const Vector & x);
+real_t diff;
+
+// We can minimize the expression |\nu \nabla u_h^* + q_h |^2 over a single element K,
+// for p+1 degree u_h^*, with the constraint \int_K u_h^* = \int_K u_h, so the mean
+// of u_h^* is the same as the one of u_h.
+//
+// This results in the problem
+//
+// (nabla w_h, \nu \nabla u_h^*) = -(nabla w_h, q_h)
+// (1, u_h^*)                    = (1, u_h)
+//
+// Since the fist equation on its own would generate a singular problem
+// the last line of the system is rewritten by the second equation.
+//
+// This elementwise operation will provide a superconvergent solution
+// \|u-u_h\|_{L^2} < C h^{p+2} |u|_{p+1}
+class HDGPostProcessing
+{
+private:
+   GridFunction *q, *u;
+
+   FiniteElementSpace *fes;
+
+   Coefficient *diffcoeff;
+
+protected:
+   const IntegrationRule *IntRule;
+
+public:
+   HDGPostProcessing(FiniteElementSpace *f, GridFunction &_q, GridFunction &_u,
+                     Coefficient &_diffcoeff)
+      : q(&_q), u(&_u), fes(f), diffcoeff(&_diffcoeff)
+   {
+      IntRule = NULL;
+   }
+
+   void Postprocessing(GridFunction &u_postprocessed) ;
+};
+
+int main(int argc, char *argv[])
+{
+   StopWatch chrono;
+
+   // 1. Parse command-line options.
+   const char *mesh_file = "../../data/inline-tri.mesh";
+   int order = 1;
+   int initial_ref_levels = 0;
+   int total_ref_levels = 2;
+   bool visualization = true;
+   bool post = true;
+   bool save = true;
+   bool hdg = true;
+   real_t memA = 0.0;
+   real_t memB = 0.0;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&initial_ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly for the initial calculation.");
+   args.AddOption(&total_ref_levels, "-tr", "--totalrefine",
+                  "Number of times to refine the mesh uniformly to get the convergence rates.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&post, "-post", "--postprocessing",
+                  "-no-post", "--no-postprocessing",
+                  "Enable or disable postprocessing.");
+   args.AddOption(&save, "-save", "--save-files", "-no-save",
+                  "--no-save-files",
+                  "Enable or disable file saving.");
+   args.AddOption(&hdg, "-hdg", "--hybrid", "-edg",
+                  "--embedded",
+                  "HDG / EDG option.");
+   args.AddOption(&memA, "-memA", "--memoryA",
+                  "Storage of A.");
+   args.AddOption(&memB, "-memB", "--memoryB",
+                  "Storage of B.");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // memA, memB \in [0,1], memB <= memA
+   if (memB > memA)
+   {
+      std::cout << "memB cannot be more than memA. Resetting to be equal" << std::endl
+                << std::flush;
+      memA = memB;
+   }
+   if (memA > 1.0)
+   {
+      std::cout << "memA cannot be more than 1. Resetting to 1" << std::endl <<
+                std::flush;
+      memA = 1.0;
+   }
+   else if (memA < 0.0)
+   {
+      std::cout << "memA cannot be less than 0. Resetting to 0." << std::endl <<
+                std::flush;
+      memA = 0.0;
+   }
+   if (memB > 1.0)
+   {
+      std::cout << "memB cannot be more than 1. Resetting to 1" << std::endl <<
+                std::flush;
+      memB = 1.0;
+   }
+   else if (memB < 0.0)
+   {
+      std::cout << "memB cannot be less than 0. Resetting to 0." << std::endl <<
+                std::flush;
+      memB = 0.0;
+   }
+
+   // 2. Read the mesh from the given mesh file. Refine it up to the initial_ref_levels.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   if (mesh->Nonconforming())
+   {
+      cout << "The current implementation does not support Nonconforming meshes. Terminating"
+           << endl << flush;
+      return 1;
+   }
+
+   for (int ii=0; ii<initial_ref_levels; ii++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 3. Vectors for the different discretization errors
+   Vector u_l2errors(total_ref_levels), q_l2errors(total_ref_levels),
+          mean_l2errors(total_ref_levels), u_star_l2errors(total_ref_levels);
+
+   // 4. Define a finite element collections and spaces on the mesh.
+   FiniteElementCollection *dg_coll(new DG_FECollection(order, dim));
+   FiniteElementCollection *face = NULL;
+   if (hdg)
+   {
+      face = new DG_Interface_FECollection(order, dim);
+   }
+   else
+   {
+      face = new H1_Trace_FECollection(order, dim);
+   }
+
+   // Finite element spaces:
+   // V_space is the vector valued DG space on elements for q_h
+   // W_space is the scalar DG space on elements for u_h
+   // M_space is the DG space on faces for lambda_h
+   FiniteElementSpace *V_space = new FiniteElementSpace(mesh, dg_coll, dim);
+   FiniteElementSpace *W_space = new FiniteElementSpace(mesh, dg_coll);
+   FiniteElementSpace *M_space = new FiniteElementSpace(mesh, face);
+
+   // 5. Define the coefficients, the exact solutions, the right hand side and the diffusion coefficient along with the diffusion penalty parameter.
+   FunctionCoefficient fcoeff(fFun);
+
+   FunctionCoefficient ucoeff(uFun_ex);
+   VectorFunctionCoefficient qcoeff(dim, qFun_ex);
+
+   diff = 1.;
+   ConstantCoefficient diffusion(diff); // diffusion constant
+   real_t tau_D = 5.0;
+
+   // 6. Define the different forms and gridfunctions.
+   HDGBilinearForm *AVarf(new HDGBilinearForm(V_space, W_space, M_space));
+   AVarf->AddHDGDomainIntegrator(new HDGDomainIntegratorDiffusion(diffusion));
+   AVarf->AddHDGFaceIntegrator(new HDGFaceIntegratorDiffusion(tau_D));
+
+   GridFunction lambda_variable(M_space);
+   GridFunction q_variable(V_space), u_variable(W_space);
+
+   LinearForm *fform(new LinearForm);
+   fform->AddDomainIntegrator(new DomainLFIntegrator(fcoeff));
+
+   for (int ref_levels = initial_ref_levels;
+        ref_levels < (initial_ref_levels + total_ref_levels); ref_levels++)
+   {
+      // 7. Compute the problem size and define the right hand side vectors
+      int dimV = V_space->GetVSize();
+      int dimW = W_space->GetVSize();
+      int dimM = M_space->GetVSize();
+
+      std::cout << "***********************************************************\n";
+      std::cout << "dim(W) = " << dimV << "\n";
+      std::cout << "dim(V) = " << dimW << "\n";
+      std::cout << "dim(M) = " << dimM << "\n";
+      std::cout << "dim(W+V+M) = " << dimV + dimW + dimM << "\n";
+      std::cout << "***********************************************************\n";
+
+      Vector rhs_R(dimV);
+      Vector rhs_F(dimW);
+      Vector V_aux(dimV);
+      Vector W_aux(dimW);
+
+      V_aux = 0.0;
+      W_aux = 0.0;
+      rhs_R = 0.0;
+
+      // 8. To eliminate the boundary conditions we project the BC to a grid function
+      // defined for the facet unknowns.
+      FunctionCoefficient lambda_coeff(uFun_ex);
+      lambda_variable.ProjectCoefficientSkeleton(lambda_coeff);
+
+      Array<int> ess_bdr(mesh->bdr_attributes.Max());
+      ess_bdr = 1;
+
+      // 9. Assemble the RHS and the Schur complement
+      fform->Update(W_space, rhs_F, 0);
+      fform->Assemble();
+
+      GridFunction *R = new GridFunction(V_space, rhs_R);
+      GridFunction *F = new GridFunction(W_space, rhs_F);
+      AVarf->AssembleSC(R, F, ess_bdr, lambda_variable, memA, memB);
+      AVarf->Finalize();
+
+      SparseMatrix* SC = AVarf->SpMatSC();
+
+      Vector* SC_RHS = AVarf->VectorSC();
+      // AVarf->VectorSC() provides -C*A^{-1} RF, the RHS for the
+      // Schur complement is  L - C*A^{-1} RF, but L is zero for this case.
+
+      // 10. Solve the Schur complement system
+      int maxIter(4000);
+      real_t rtol(1.e-13);
+      real_t atol(0.0);
+      GSSmoother M(*SC);
+      BiCGSTABSolver solver;
+      solver.SetAbsTol(atol);
+      solver.SetRelTol(rtol);
+      solver.SetMaxIter(maxIter);
+      solver.SetOperator(*SC);
+      solver.SetPrintLevel(-1);
+      solver.SetPreconditioner(M);
+      chrono.Clear();
+      chrono.Start();
+      solver.Mult(*SC_RHS, lambda_variable);
+      chrono.Stop();
+
+      if (solver.GetConverged())
+         std::cout << "Iterative method converged in " << solver.GetNumIterations()
+                   << " iterations with a residual norm of " << solver.GetFinalNorm() << ".\n";
+      else
+         std::cout << "Iterative method did not converge in " <<
+                   solver.GetNumIterations()
+                   << " iterations. Residual norm is " << solver.GetFinalNorm() << ".\n";
+      std::cout << "Iterative method solver took " << chrono.RealTime() << "s. \n";
+
+      // Delete the SC matrix to save memory
+      SC = NULL;
+
+      // 11. Reconstruction
+      // Reconstruct the solution u and q from the facet solution lambda
+      AVarf->Reconstruct(R, F, &lambda_variable, &q_variable, &u_variable);
+
+      // 12. Compute the discretization error
+      int order_quad = max(2, 2*order+2);
+      const IntegrationRule *irs[Geometry::NumGeom];
+      for (int i=0; i < Geometry::NumGeom; ++i)
+      {
+         irs[i] = &(IntRules.Get(i, order_quad));
+      }
+      real_t err_u    = u_variable.ComputeL2Error(ucoeff, irs);
+      real_t err_q    = q_variable.ComputeL2Error(qcoeff, irs);
+      real_t err_mean  = u_variable.ComputeMeanLpError(2.0, ucoeff, irs);
+
+      u_l2errors(ref_levels-initial_ref_levels) = fabs(err_u);
+      q_l2errors(ref_levels-initial_ref_levels) = fabs(err_q);
+      mean_l2errors(ref_levels-initial_ref_levels) = fabs(err_mean);
+
+      std::cout << "|| u_h - u_ex || = " << err_u << "\n";
+      std::cout << "|| q_h - q_ex || = " << err_q << "\n";
+      std::cout << "|| mean(u_h) - mean(u_ex) || = " << err_mean << "\n";
+
+      // 13. Save the mesh and the solution.
+      if (save)
+      {
+         ofstream mesh_ofs("ex_hdg.mesh");
+         mesh_ofs.precision(8);
+         mesh->Print(mesh_ofs);
+
+         ofstream q_variable_ofs("sol_q.gf");
+         q_variable_ofs.precision(8);
+         q_variable.Save(q_variable_ofs);
+
+         ofstream u_variable_ofs("sol_u.gf");
+         u_variable_ofs.precision(8);
+         u_variable.Save(u_variable_ofs);
+
+         ofstream lambda_variable_ofs("sol_lambda.gf");
+         lambda_variable_ofs.precision(8);
+         lambda_variable.Save(lambda_variable_ofs);
+      }
+
+      // 14. Send the solution by socket to a GLVis server.
+      if (visualization)
+      {
+         char vishost[] = "localhost";
+         int  visport   = 19916;
+         socketstream u_sock(vishost, visport);
+         u_sock.precision(8);
+         u_sock << "solution\n" << *mesh << u_variable << "window_title 'Solution u'" <<
+                endl;
+
+         socketstream q_sock(vishost, visport);
+         q_sock.precision(8);
+         q_sock << "solution\n" << *mesh << q_variable << "window_title 'Solution q'" <<
+                endl;
+      }
+
+      // 15. Postprocessing
+      if (post)
+      {
+         FiniteElementCollection *dg_coll_pstar(new DG_FECollection(order+1, dim));
+         FiniteElementSpace *Vstar_space = new FiniteElementSpace(mesh, dg_coll_pstar);
+
+         GridFunction u_post(Vstar_space);
+
+         HDGPostProcessing *hdgpost(new HDGPostProcessing(Vstar_space, q_variable,
+                                                          u_variable, diffusion));
+
+         hdgpost->Postprocessing(u_post);
+
+         order_quad = max(2, 2*order+5);
+         for (int i=0; i < Geometry::NumGeom; ++i)
+         {
+            irs[i] = &(IntRules.Get(i, order_quad));
+         }
+         real_t err_u_post   = u_post.ComputeL2Error(ucoeff, irs);
+
+         u_star_l2errors(ref_levels-initial_ref_levels) = fabs(err_u_post);
+
+         std::cout << "|| u^*_h - u_ex || = " << err_u_post << "\n";
+
+         if (save)
+         {
+            ofstream u_post_ofs("sol_u_star.gf");
+            u_post_ofs.precision(8);
+            u_post.Save(u_post_ofs);
+         }
+
+         if (visualization)
+         {
+            char vishost[] = "localhost";
+            int  visport   = 19916;
+            socketstream u_star_sock(vishost, visport);
+            u_star_sock.precision(8);
+            u_star_sock << "solution\n" << *mesh << u_post <<
+                        "window_title 'Solution u_star'" << endl;
+         }
+
+         delete hdgpost;
+         delete Vstar_space;
+         delete dg_coll_pstar;
+      }
+
+      // 16. Refine the mesh to increase the resolution and update the spaces and the forms.
+      mesh->UniformRefinement();
+
+      V_space->Update(0);
+      W_space->Update(0);
+      M_space->Update(0);
+
+      AVarf->Update();
+      q_variable.Update();
+      u_variable.Update();
+      lambda_variable.Update();
+
+      delete R;
+      delete F;
+   }
+
+   // 17. Print the results
+   std::cout << "\n\n-----------------------\n";
+   std::cout <<
+             "level  u_l2errors  order   q_l2errors  order   mean_l2errors  order u_star_l2errors   order\n";
+   std::cout << "-----------------------\n";
+   for (int ref_levels = 0; ref_levels < total_ref_levels; ref_levels++)
+   {
+      if (ref_levels == 0)
+      {
+         std::cout << "  " << ref_levels << "    "
+                   << std::setprecision(2) << std::scientific << u_l2errors(ref_levels)
+                   << "   " << " -       "
+                   << std::setprecision(2) << std::scientific << q_l2errors(ref_levels)
+                   << "    " << " -       "
+                   << std::setprecision(2) << std::scientific << mean_l2errors(ref_levels)
+                   << "    " << " -       "
+                   << std::setprecision(2) << std::scientific << u_star_l2errors(ref_levels)
+                   << "    " << " -       " << std::endl;
+      }
+      else
+      {
+         real_t u_order     = log(u_l2errors(ref_levels)/u_l2errors(ref_levels-1))/log(
+                                 0.5);
+         real_t q_order     = log(q_l2errors(ref_levels)/q_l2errors(ref_levels-1))/log(
+                                 0.5);
+         real_t mean_order   = log(mean_l2errors(ref_levels)/mean_l2errors(
+                                      ref_levels-1))/log(0.5);
+         real_t u_star_order = log(u_star_l2errors(ref_levels)/u_star_l2errors(
+                                      ref_levels-1))/log(0.5);
+         std::cout << "  " << ref_levels << "    "
+                   << std::setprecision(2) << std::scientific << u_l2errors(ref_levels)
+                   << "  " << std::setprecision(4) << std::fixed << u_order
+                   << "    " << std::setprecision(2) << std::scientific << q_l2errors(ref_levels)
+                   << "   " << std::setprecision(4) << std::fixed << q_order
+                   << "    " << std::setprecision(2) << std::scientific << mean_l2errors(
+                      ref_levels)
+                   << "   " << std::setprecision(4) << std::fixed << mean_order
+                   << "    " << std::setprecision(2) << std::scientific << u_star_l2errors(
+                      ref_levels)
+                   << "   " << std::setprecision(4) << std::fixed << u_star_order << std::endl;
+      }
+   }
+   std::cout << "\n\n";
+
+   // 18. Free the used memory.
+   delete mesh;
+   delete V_space;
+   delete W_space;
+   delete M_space;
+   delete AVarf;
+   delete fform;
+   delete dg_coll;
+   delete face;
+
+   std::cout << "Done." << std::endl ;
+
+   return 0;
+}
+
+
+real_t uFun_ex(const Vector & x)
+{
+   real_t xi(x(0));
+   real_t yi(x(1));
+
+   int dim = x.Size();
+
+   switch (dim)
+   {
+      case 2:
+      {
+         return 1.0 + xi + sin(2.0*M_PI*xi)*sin(2.0*M_PI*yi);
+         break;
+      }
+      case 3:
+      {
+         real_t zi(x(2));
+         return xi + sin(2.0*M_PI*xi)*sin(2.0*M_PI*yi)*sin(2.0*M_PI*zi);
+         break;
+      }
+   }
+
+   return 0;
+}
+
+void qFun_ex(const Vector & x, Vector & q)
+{
+   real_t xi(x(0));
+   real_t yi(x(1));
+   int dim = x.Size();
+
+   switch (dim)
+   {
+      case 2:
+      {
+         q(0) = -diff*1.0 - diff*2.0*M_PI*cos(2.0*M_PI*xi)*sin(2.0*M_PI*yi);
+         q(1) =  0.0 - diff*2.0*M_PI*sin(2.0*M_PI*xi)*cos(2.0*M_PI*yi);
+         break;
+      }
+      case 3:
+      {
+         real_t zi(x(2));
+         q(0) = -diff*1.0 - diff*2.0*M_PI*cos(2.0*M_PI*xi)*sin(2.0*M_PI*yi)*sin(
+                   2.0*M_PI*zi);
+         q(1) =  0.0 - diff*2.0*M_PI*sin(2.0*M_PI*xi)*cos(2.0*M_PI*yi)*sin(2.0*M_PI*zi);
+         q(2) =  0.0 - diff*2.0*M_PI*sin(2.0*M_PI*xi)*sin(2.0*M_PI*yi)*cos(2.0*M_PI*zi);
+         break;
+      }
+   }
+}
+
+
+real_t fFun(const Vector & x)
+{
+   real_t xi(x(0));
+   real_t yi(x(1));
+   int dim = x.Size();
+
+   switch (dim)
+   {
+      case 2:
+      {
+         return diff*8.0*M_PI*M_PI*sin(2.0*M_PI*xi)*sin(2.0*M_PI*yi);
+         break;
+      }
+      case 3:
+      {
+         real_t zi(x(2));
+         return diff*12.0*M_PI*M_PI*sin(2.0*M_PI*xi)*sin(2.0*M_PI*yi)*sin(2.0*M_PI*zi);
+         break;
+      }
+   }
+
+   return 0;
+
+}
+
+// Postprocessing
+void HDGPostProcessing::Postprocessing(GridFunction &u_postprocessed)
+{
+   Mesh *mesh = fes->GetMesh();
+   Array<int>  vdofs;
+   Vector      elmat2, shape, RHS, to_RHS, vals, uvals;
+   real_t      RHS2;
+   DenseMatrix elmat, invdfdx, dshape, dshapedxt, qvals;
+
+   int  ndofs;
+   const FiniteElement *fe_elem;
+   ElementTransformation *Trans;
+
+   for (int i = 0; i < fes->GetNE(); i++)
+   {
+      fes->GetElementVDofs(i, vdofs);
+      ndofs = vdofs.Size();
+      vals.SetSize(ndofs);
+      // elmat is the matrix for the -(nabla w_h, q_h) term
+      elmat.SetSize(ndofs);
+      // elmat 1 is the vector for the (1, u_h^*) term
+      elmat2.SetSize(ndofs);
+      shape.SetSize(ndofs);
+
+      RHS.SetSize(ndofs);
+      to_RHS.SetSize(ndofs);
+
+      elmat = 0.0;
+      elmat2 = 0.0;
+      RHS = 0.0;
+      RHS2 = 0.0;
+
+      fe_elem = fes->GetFE(i);
+      int dim = fe_elem->GetDim();
+      int spaceDim = dim;
+      invdfdx.SetSize(dim, spaceDim);
+      dshape.SetSize(ndofs, spaceDim);
+      dshapedxt.SetSize(ndofs, spaceDim);
+
+      Trans = mesh->GetElementTransformation(i);
+
+      const IntegrationRule *ir = IntRule;
+      if (ir == NULL)
+      {
+         int order = 3*fe_elem->GetOrder() + 3;
+         ir = &IntRules.Get(fe_elem->GetGeomType(), order);
+      }
+
+      // Get the values of u_h and q_h
+      u->GetValues(i, *ir, uvals);
+      q->GetVectorValues(*Trans, *ir, qvals);
+
+      for (int j = 0; j < ir->GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(j);
+
+         fe_elem->CalcDShape(ip, dshape);
+         fe_elem->CalcShape(ip, shape);
+
+         Trans->SetIntPoint(&ip);
+         // Compute invdfdx = / adj(J),       if J is square
+         //               \ adj(J^t.J).J^t, otherwise
+         CalcAdjugate(Trans->Jacobian(), invdfdx);
+         real_t w = Trans->Weight();
+         w = ip.weight / w;
+         w *= diffcoeff->Eval(*Trans, ip);
+         Mult(dshape, invdfdx, dshapedxt);
+
+         // compute the (nabla w_h, \nu \nabla u_h^*) term
+         AddMult_a_AAt(w, dshapedxt, elmat);
+
+         dshapedxt *= ip.weight ;
+
+         Vector qval_col;
+         qvals.GetColumn(j, qval_col);
+
+         // compute (nabla w_h, q_h)
+         dshapedxt.Mult(qval_col, to_RHS);
+
+         // subtract it from the rhs
+         RHS -= to_RHS;
+
+         // compute (1, u_h^*)
+         shape *= (Trans->Weight() * ip.weight);
+         elmat2 += shape;
+
+         // compute (1, u_h)
+         real_t rhs_weight = (Trans->Weight() * ip.weight);
+         RHS2  += (uvals(j)*rhs_weight);
+
+      }
+
+      // changing the last row and the last entry
+      for (int j = 0; j < ndofs; j++)
+      {
+         elmat(ndofs-1,j) = elmat2(j);
+      }
+      RHS(ndofs-1) = RHS2;
+
+      // solve the local problem
+      elmat.Invert();
+      elmat.Mult(RHS, vals);
+      u_postprocessed.SetSubVector(vdofs, vals);
+
+   }
+}

miniapps/hdg/hdg_poissonp.cpp

@@ -0,0 +1,831 @@
+//                  MFEM Example Hybridizable DG
+//
+// Compile with: make hdg_poissonp
+//
+// Sample runs:  mpirun -np 1 hdg_poissonp -o 1 -r 1 -tr 4 -no-vis
+//         mpirun -np 2 hdg_poissonp -o 5 -r 1 -tr 4 -no-vis
+//         mpirun -np 2 hdg_poissonp -o 1 -r 4 -tr 1
+//         mpirun -np 3 hdg_poissonp -o 5 -r 4 -tr 1
+//         mpirun -np 2 hdg_poissonp -o 1 -r 1 -tr 4 -no-vis -m ../data/inline-tri.mesh
+//         mpirun -np 2 hdg_poissonp -o 5 -r 1 -tr 4 -no-vis -m ../data/inline-tri.mesh
+//         mpirun -np 4 hdg_poissonp -o 1 -r 5 -tr 1 -m ../data/inline-tri.mesh
+//         mpirun -np 2 hdg_poissonp -o 5 -r 5 -tr 1 -m ../data/inline-tri.mesh
+//
+// Description:  This example code solves the 2D/3D diffusion problem
+//                -\nu Delta u = f
+//               with Dirichlet boundary conditions, using HDG discretization.
+//
+// The methods approximates the solution u, the diffusive flux q = -\nu \nabla u,
+// and the restriction of u to the faces, denoted by lambda.
+//
+// The weak form is: seek (q,u,\lambda) such that for all (v, w, \mu)
+//
+// -\nu^{-1}(q, v)       + (u, div(v))       - <\lambda, v \cdot n>      = 0
+//  (div(q), w)          + <\tau u, w>       - <\tau \lambda, w>         = (f, w)
+// -<[[q \cdot n]], \mu> - <[[\tau u]], \mu> + <[[(\tau \lambda]], \mu>  = 0
+//
+// where [[.]] is the jump operator, (.,.) is the d-dimensional L2 product,
+// <.,.> is the d-1 dimensional L2 product.
+//
+// The discretization is based on the paper:
+//
+// N.C. Nguyen, J. Peraire, B. Cockburn, An implicit high-order hybridizable
+// discontinuous Galerkin method for linear convection–diffusion equations,
+// J. Comput. Phys., 2009, 228:9, 3232--3254.
+//
+// Contributed by: T. Horvath, S. Rhebergen, A. Sivas
+//                 University of Waterloo
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include <algorithm>
+#include "HDGBilinearForm.hpp"
+#include "hdg_integrators.hpp"
+
+using namespace std;
+using namespace mfem;
+
+// Define the analytical solution and forcing terms / boundary conditions
+real_t uFun_ex(const Vector & x);
+void qFun_ex(const Vector & x, Vector & q);
+real_t fFun(const Vector & x);
+real_t diff;
+
+// We can minimize the expression |\nu \nabla u_h^* + q_h |^2 over a single element K,
+// for p+1 degree u_h^*, with the constraint \int_K u_h^* = \int_K u_h, so the mean
+// of u_h^* is the same as the one of u_h.
+//
+// This results in the problem
+//
+// (nabla w_h, \nu \nabla u_h^*) = -(nabla w_h, q_h)
+// (1, u_h^*)                    = (1, u_h)
+//
+// Since the fist equation on its own would generate a singular problem
+// the last line of the system is rewritten by the second equation.
+//
+// This elementwise operation will provide a superconvergent solution
+// \|u-u_h\|_{L^2} < C h^{p+2} |u|_{p+1}
+class pHDGPostProcessing
+{
+private:
+   ParGridFunction *q, *u;
+
+   ParFiniteElementSpace *pfes;
+
+   Coefficient *diffcoeff;
+protected:
+   const IntegrationRule *IntRule;
+
+public:
+   pHDGPostProcessing(ParFiniteElementSpace *f, ParGridFunction &_q,
+                      ParGridFunction &_u, Coefficient &_diffcoeff)
+      : q(&_q), u(&_u), pfes(f), diffcoeff(&_diffcoeff)
+   {
+      IntRule = NULL;
+   }
+
+   void Postprocessing(ParGridFunction &u_postprocessed) ;
+};
+
+int main(int argc, char *argv[])
+{
+   StopWatch chrono;
+
+   // 1. Initialize MPI.
+   Mpi::Init(argc, argv);
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();
+
+   real_t assemblyTime, solveTime, reconstructTime, pprocessTime;
+   real_t GassemblyTime, GsolveTime, GreconstructTime, GpprocessTime;
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../../data/inline-tri.mesh";
+   int order = 1;
+   int initial_ref_levels = 0;
+   int total_ref_levels = 2;
+   bool visualization = true;
+   bool verbose = (myid == 0);
+   bool post = true;
+   bool save = true;
+   bool hdg = true;
+   real_t memA = 0.0;
+   real_t memB = 0.0;
+   bool petsc = false;
+   const char *petscrc_file = "";
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&initial_ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly for the initial calculation.");
+   args.AddOption(&total_ref_levels, "-tr", "--totalrefine",
+                  "Number of times to refine the mesh uniformly to get the convergence rates.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&post, "-post", "--postprocessing",
+                  "-no-post", "--no-postprocessing",
+                  "Enable or disable postprocessing.");
+   args.AddOption(&save, "-save", "--save-files", "-no-save",
+                  "--no-save-files",
+                  "Enable or disable file saving.");
+   args.AddOption(&hdg, "-hdg", "--hybrid", "-edg",
+                  "--embedded",
+                  "HDG / EDG option.");
+   args.AddOption(&memA, "-memA", "--memoryA",
+                  "Storage of A.");
+   args.AddOption(&memB, "-memB", "--memoryB",
+                  "Storage of B.");
+   args.AddOption(&petsc, "-petsc", "--use-petsc",
+                  "-no-petsc", "--no-use-petsc",
+                  "Enable or disable SC solver.");
+   args.AddOption(&petscrc_file, "-petscopts", "--petscopts",
+                  "PetscOptions file to use.");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      if (verbose)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (verbose)
+   {
+      args.PrintOptions(cout);
+   }
+
+#ifdef MFEM_USE_PETSC
+   // We initialize PETSc
+   MFEMInitializePetsc(NULL,NULL,petscrc_file,NULL);
+#endif
+
+#ifndef MFEM_USE_PETSC
+   if (petsc)
+   {
+      std::cout << "MFEM does not use PETSc. Change the solver to hypre" << std::endl
+                << std::flush;
+      petsc = false;
+   }
+#endif
+
+   // memA, memB \in [0,1], memB <= memA
+   if (memB > memA)
+   {
+      std::cout << "memB cannot be more than memA. Resetting to be equal" << std::endl
+                << std::flush;
+      memA = memB;
+   }
+   if (memA > 1.0)
+   {
+      std::cout << "memA cannot be more than 1. Resetting to 1" << std::endl <<
+                std::flush;
+      memA = 1.0;
+   }
+   else if (memA < 0.0)
+   {
+      std::cout << "memA cannot be less than 0. Resetting to 0." << std::endl <<
+                std::flush;
+      memA = 0.0;
+   }
+   if (memB > 1.0)
+   {
+      std::cout << "memB cannot be more than 1. Resetting to 1" << std::endl <<
+                std::flush;
+      memB = 1.0;
+   }
+   else if (memB < 0.0)
+   {
+      std::cout << "memB cannot be less than 0. Resetting to 0." << std::endl <<
+                std::flush;
+      memB = 0.0;
+   }
+
+   // 3. Read the mesh from the given mesh file. Refine it up to the initial_ref_levels.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   if (mesh->Nonconforming())
+   {
+      if (verbose)
+      {
+         cout << "The current implementation does not support Nonconforming meshes. Terminating"
+              << endl << flush;
+      }
+#ifdef MFEM_USE_PETSC
+      MFEMFinalizePetsc();
+#endif
+      MPI_Finalize();
+      return 1;
+   }
+
+   for (int ii=0; ii<initial_ref_levels; ii++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+
+   // 4. Vectors for the different discretization errors
+   Vector u_l2errors(total_ref_levels), q_l2errors(total_ref_levels),
+          mean_l2errors(total_ref_levels), u_star_l2errors(total_ref_levels);
+
+   // 5. Define a finite element collections and spaces on the mesh.
+   FiniteElementCollection *dg_coll(new DG_FECollection(order, dim));
+   FiniteElementCollection *face = NULL;
+   if (hdg)
+   {
+      face = new DG_Interface_FECollection(order, dim);
+   }
+   else
+   {
+      face = new H1_Trace_FECollection(order, dim);
+   }
+
+   // Finite element spaces:
+   // V_space is the vector valued DG space on elements for q_h
+   // W_space is the scalar DG space on elements for u_h
+   // M_space is the DG space on faces for lambda_h
+   ParFiniteElementSpace *V_space = new ParFiniteElementSpace(pmesh, dg_coll, dim);
+   ParFiniteElementSpace *W_space = new ParFiniteElementSpace(pmesh, dg_coll);
+   ParFiniteElementSpace *M_space = new ParFiniteElementSpace(pmesh, face);
+
+   // 6. Define the coefficients, the exact solutions, the right hand side and the diffusion coefficient along with the diffusion penalty parameter.
+   FunctionCoefficient fcoeff(fFun);
+
+   FunctionCoefficient ucoeff(uFun_ex);
+   VectorFunctionCoefficient qcoeff(dim, qFun_ex);
+
+   diff = 1.;
+   ConstantCoefficient diffusion(diff); // diffusion constant
+   real_t tau_D = 5.0;
+
+   // 7. Define the different forms and gridfunctions.
+   HDGBilinearForm *AVarf(new HDGBilinearForm(V_space, W_space, M_space, true));
+   AVarf->AddHDGDomainIntegrator(new HDGDomainIntegratorDiffusion(diffusion));
+   AVarf->AddHDGFaceIntegrator(new HDGFaceIntegratorDiffusion(tau_D));
+
+   ParGridFunction lambda(M_space);
+   ParGridFunction q_variable(V_space), u_variable(W_space);
+
+   ParLinearForm *fform(new ParLinearForm);
+   fform->AddDomainIntegrator(new DomainLFIntegrator(fcoeff));
+
+   for (int ref_levels = initial_ref_levels;
+        ref_levels < (initial_ref_levels + total_ref_levels); ref_levels++)
+   {
+      // 8. Compute the problem size and define the right hand side vectors
+      HYPRE_Int dimV = V_space->GlobalTrueVSize();
+      HYPRE_Int dimW = W_space->GlobalTrueVSize();
+      HYPRE_Int dimM = M_space->GlobalTrueVSize();
+
+      if (verbose)
+      {
+         std::cout << "***********************************************************\n";
+         std::cout << "dim(V) = " << dimV << "\n";
+         std::cout << "dim(W) = " << dimW << "\n";
+         std::cout << "dim(M) = " << dimM << "\n";
+         std::cout << "dim(V+W+M) = " << dimV + dimW + dimM << "\n";
+         std::cout << "***********************************************************\n";
+      }
+
+      HypreParVector *trueR(new HypreParVector(V_space));
+      *trueR = 0.0;
+
+      HypreParVector *trueF;
+
+      // 9. To eliminate the boundary conditions we project the BC to a grid function
+      // defined for the facet unknowns.
+      lambda.ProjectCoefficientSkeleton(ucoeff);
+
+      HypreParVector *Lambda =  new HypreParVector(M_space);
+      lambda.ParallelProject(*Lambda);
+
+      Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+      ess_bdr = 1;
+
+      // 10. Assemble the RHS and the Schur complement
+      fform->Update(W_space);
+      fform->Assemble();
+
+      trueF = fform->ParallelAssemble();
+
+      // Creating a gridfunctions for the elimination of the boundary
+      ParGridFunction *R = new ParGridFunction(V_space, trueR);
+      ParGridFunction *F = new ParGridFunction(W_space, trueF);
+
+      chrono.Clear();
+      chrono.Start();
+      AVarf->AssembleSC(R, F, ess_bdr, lambda, memA, memB);
+      chrono.Stop();
+      AVarf->Finalize();
+
+      assemblyTime = chrono.RealTime();
+      HypreParMatrix *SC = AVarf->ParallelAssembleSC();
+
+      HypreParVector *rhs_SC = AVarf->ParallelVectorSC();
+      // AVarf->ParallelVectorSC() provides -C*A^{-1} RF, the RHS for the
+      // Schur complement is  L - C*A^{-1} RF, but L is zero for this case
+
+      // 11. Solve the Schur complement system
+      real_t tol = 1.0e-12;
+      int maxIter = 1000;
+      int PrintLevel = -1;
+
+      if (petsc)
+      {
+#ifdef MFEM_USE_PETSC
+         // Solver using PETSc
+         //=======================
+         PetscLinearSolver *petsc_solver;
+         PetscPreconditioner *petsc_precon= NULL;
+         petsc_solver = new PetscLinearSolver(MPI_COMM_WORLD, "solver_", 1, 0);
+         petsc_precon = new PetscPreconditioner(MPI_COMM_WORLD,*SC,"solver_");
+         petsc_solver->SetOperator(*SC);
+         petsc_solver->SetPreconditioner(*petsc_precon);
+         petsc_solver->SetTol(tol);
+         petsc_solver->SetAbsTol(0.0);
+         petsc_solver->SetMaxIter(maxIter);
+         petsc_solver->SetPrintLevel(PrintLevel);
+         chrono.Clear();
+         chrono.Start();
+         petsc_solver->Mult(*rhs_SC, *Lambda);
+         chrono.Stop();
+
+         if (verbose)
+         {
+            if (petsc_solver->GetConverged())
+               std::cout << "Solver converged in " << petsc_solver->GetNumIterations()
+                         << " iterations with a residual norm of " << petsc_solver->GetFinalNorm() <<
+                         ".\n";
+            else
+               std::cout << "Solver did not converge in " << petsc_solver->GetNumIterations()
+                         << " iterations. Residual norm is " << petsc_solver->GetFinalNorm() << ".\n";
+            std::cout << "Solver solver took " << chrono.RealTime() << "s. \n";
+         }
+
+         delete petsc_solver;
+         delete petsc_precon;
+#endif
+      }
+      else
+      {
+         HypreBoomerAMG *amg = new HypreBoomerAMG(*SC);
+         HyprePCG *pcg = new HyprePCG(*SC);
+         pcg->SetTol(tol);
+         pcg->SetMaxIter(maxIter);
+         amg->SetPrintLevel(PrintLevel);
+         pcg->SetPrintLevel(PrintLevel);
+         pcg->SetPreconditioner(*amg);
+         chrono.Clear();
+         chrono.Start();
+         pcg->Mult(*rhs_SC, *Lambda);
+         chrono.Stop();
+
+         int numIterations = 0;
+         pcg->GetNumIterations(numIterations);
+
+         if (verbose)
+         {
+            std::cout << "\nIterative method converged in "
+                      << numIterations << ".\n";
+
+            std::cout << "Iterative solver took " << chrono.RealTime() << "s. \n";
+         }
+      }
+
+      // Delete the SC matrix to save memory
+      SC = NULL;
+      solveTime = chrono.RealTime();
+
+      // 12. Reconstruction
+      // Create a gridfunction from the right hand side.
+      // It is mostly important for the parallel code,
+      // here it is done this way to make the 2 codes more similar
+      lambda = ParGridFunction(M_space, Lambda);
+
+      chrono.Clear();
+      chrono.Start();
+      AVarf->Reconstruct(R, F, &lambda, &q_variable, &u_variable);
+      chrono.Stop();
+
+      reconstructTime = chrono.RealTime();
+
+      // 13. Compute the discretization error
+      int order_quad = max(2, 2*order+1);
+      const IntegrationRule *irs[Geometry::NumGeom];
+      for (int i=0; i < Geometry::NumGeom; ++i)
+      {
+         irs[i] = &(IntRules.Get(i, order_quad));
+      }
+
+      real_t err_u   = u_variable.ComputeL2Error(ucoeff, irs);
+      real_t norm_u   = ComputeGlobalLpNorm(2., ucoeff, *pmesh, irs);
+      real_t err_q   = q_variable.ComputeL2Error(qcoeff, irs);
+      real_t norm_q   = ComputeGlobalLpNorm(2., qcoeff, *pmesh, irs);
+      real_t err_mean  = u_variable.ComputeMeanLpError(2.0, ucoeff, irs);
+
+      if (verbose)
+      {
+         std::cout << "|| u_h - u_ex || / || u_ex || = " << err_u / norm_u << "\n";
+         std::cout << "|| q_h - q_ex || / || q_ex || = " << err_q / norm_q << "\n";
+         std::cout << "|| u_h - u_ex || = " << err_u << "\n";
+         std::cout << "|| q_h - q_ex || = " << err_q << "\n";
+         std::cout << "|| mean(u_h) - mean(u_ex) || = " << err_mean << "\n";
+      }
+
+      u_l2errors(ref_levels-initial_ref_levels) = fabs(err_u);
+      q_l2errors(ref_levels-initial_ref_levels) = fabs(err_q);
+      mean_l2errors(ref_levels-initial_ref_levels) = fabs(err_mean);
+
+
+      // 14. Save the mesh and the solution.
+      if (save)
+      {
+         ostringstream mesh_name, u_name, q_name, lambda_name;
+         mesh_name << "mesh." << setfill('0') << setw(6) << myid;
+         u_name << "sol_u." << setfill('0') << setw(6) << myid;
+         q_name << "sol_q." << setfill('0') << setw(6) << myid;
+         lambda_name << "sol_lambda." << setfill('0') << setw(6) << myid;
+
+         ofstream mesh_ofs(mesh_name.str().c_str());
+         mesh_ofs.precision(8);
+         pmesh->Print(mesh_ofs);
+
+         ofstream u_ofs(u_name.str().c_str());
+         u_ofs.precision(8);
+         u_variable.Save(u_ofs);
+
+         ofstream q_ofs(q_name.str().c_str());
+         q_ofs.precision(8);
+         q_variable.Save(q_ofs);
+
+         ParGridFunction lambda_variable(M_space, Lambda);
+         ofstream lambda_ofs(lambda_name.str().c_str());
+         lambda_ofs.precision(8);
+         lambda_variable.Save(lambda_ofs);
+      }
+
+      // 15. Send the solution by socket to a GLVis server.
+      if (visualization)
+      {
+         char vishost[] = "localhost";
+         int  visport   = 19916;
+         socketstream u_sock(vishost, visport);
+         u_sock << "parallel " << num_procs << " " << myid << "\n";
+         u_sock.precision(8);
+         u_sock << "solution\n" << *pmesh << u_variable << "window_title 'U'"
+                << endl;
+         // Make sure all ranks have sent their 'u' solution before initiating
+         // another set of GLVis connections (one from each rank):
+         MPI_Barrier(pmesh->GetComm());
+         socketstream q_sock(vishost, visport);
+         q_sock << "parallel " << num_procs << " " << myid << "\n";
+         q_sock.precision(8);
+         q_sock << "solution\n" << *pmesh << q_variable << "window_title 'Q'"
+                << endl;
+      }
+
+      // 16. Postprocessing
+      if (post)
+      {
+         FiniteElementCollection *dg_coll_pstar(new DG_FECollection(order+1, dim));
+         ParFiniteElementSpace *Vstar_space = new ParFiniteElementSpace(pmesh,
+                                                                        dg_coll_pstar);
+
+         ParGridFunction u_post(Vstar_space);
+
+         pHDGPostProcessing *hdgpost(new pHDGPostProcessing(Vstar_space, q_variable,
+                                                            u_variable, diffusion));
+
+         chrono.Clear();
+         chrono.Start();
+         hdgpost->Postprocessing(u_post);
+         chrono.Stop();
+
+         pprocessTime = chrono.RealTime();
+
+         order_quad = max(2, 2*order+5);
+         for (int i=0; i < Geometry::NumGeom; ++i)
+         {
+            irs[i] = &(IntRules.Get(i, order_quad));
+         }
+         real_t err_u_post   = u_post.ComputeL2Error(ucoeff, irs);
+
+         u_star_l2errors(ref_levels-initial_ref_levels) = fabs(err_u_post);
+
+         if (verbose)
+         {
+            std::cout << "|| u^*_h - u_ex || = " << err_u_post << "\n";
+         }
+
+         if (save)
+         {
+            ostringstream u_star_name;
+            u_star_name << "sol_u_star." << setfill('0') << setw(6) << myid;
+            ofstream u_star_ofs(u_star_name.str().c_str());
+            u_star_ofs.precision(8);
+            u_post.Save(u_star_ofs);
+         }
+
+         if (visualization)
+         {
+            char vishost[] = "localhost";
+            int  visport   = 19916;
+
+            MPI_Barrier(pmesh->GetComm());
+            socketstream u_star_sock(vishost, visport);
+            u_star_sock << "parallel " << num_procs << " " << myid << "\n";
+            u_star_sock.precision(8);
+            u_star_sock << "solution\n" << *pmesh << u_post << "window_title 'U_star'"
+                        << endl;
+         }
+
+         delete hdgpost;
+         delete Vstar_space;
+         delete dg_coll_pstar;
+      }
+
+      // 17. Refine the mesh to increase the resolution and update the spaces and the forms. Print the runtimes
+      pmesh->UniformRefinement();
+
+      V_space->Update(0);
+      W_space->Update(0);
+      M_space->Update(0);
+
+      AVarf->Update();
+      q_variable.Update();
+      u_variable.Update();
+      lambda.Update();
+      MPI_Reduce(&assemblyTime,&GassemblyTime,1,MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD);
+      MPI_Reduce(&solveTime,&GsolveTime,1,MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD);
+      MPI_Reduce(&reconstructTime,&GreconstructTime,1,MPI_DOUBLE,MPI_MAX,0,
+                 MPI_COMM_WORLD);
+      MPI_Reduce(&pprocessTime,&GpprocessTime,1,MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD);
+
+      if (verbose)
+      {
+         printf("\t Assembly time   = %.2f\n",GassemblyTime);
+         printf("\t Solve time     = %.2f\n",GsolveTime);
+         printf("\t Reconstruct time = %.2f\n",GreconstructTime);
+         printf("\t Postprocess time = %.2f\n",GpprocessTime);
+      }
+
+      delete R;
+      delete F;
+   }
+
+   // 18. Print the results
+   if (verbose)
+   {
+      std::cout << "\n\n-----------------------\n";
+      std::cout <<
+                "level  u_l2errors  order   q_l2errors  order   mean_l2errors  order u_star_l2errors   order\n";
+      std::cout << "-----------------------\n";
+      for (int ref_levels = 0; ref_levels < total_ref_levels; ref_levels++)
+      {
+         if (ref_levels == 0)
+         {
+            std::cout << "  " << ref_levels << "   "
+                      << std::setprecision(2) << std::scientific << u_l2errors(ref_levels)
+                      << "   " << " -      "
+                      << std::setprecision(2) << std::scientific << q_l2errors(ref_levels)
+                      << "    " << " -      "
+                      << std::setprecision(2) << std::scientific << mean_l2errors(ref_levels)
+                      << "    " << " -      "
+                      << std::setprecision(2) << std::scientific << u_star_l2errors(ref_levels)
+                      << "    " << " -      " << std::endl;
+         }
+         else
+         {
+            real_t u_order    = log(u_l2errors(ref_levels)/u_l2errors(ref_levels-1))/log(
+                                   0.5);
+            real_t q_order    = log(q_l2errors(ref_levels)/q_l2errors(ref_levels-1))/log(
+                                   0.5);
+            real_t mean_order   = log(mean_l2errors(ref_levels)/mean_l2errors(
+                                         ref_levels-1))/log(0.5);
+            real_t u_star_order = log(u_star_l2errors(ref_levels)/u_star_l2errors(
+                                         ref_levels-1))/log(0.5);
+            std::cout << "  " << ref_levels << "   "
+                      << std::setprecision(2) << std::scientific << u_l2errors(ref_levels)
+                      << "  " << std::setprecision(4) << std::fixed << u_order
+                      << "   " << std::setprecision(2) << std::scientific << q_l2errors(ref_levels)
+                      << "   " << std::setprecision(4) << std::fixed << q_order
+                      << "   " << std::setprecision(2) << std::scientific << mean_l2errors(ref_levels)
+                      << "   " << std::setprecision(4) << std::fixed << mean_order
+                      << "   " << std::setprecision(2) << std::scientific << u_star_l2errors(
+                         ref_levels)
+                      << "   " << std::setprecision(4) << std::fixed << u_star_order << std::endl;
+         }
+      }
+   }
+
+   // 19. Free the used memory.
+   delete pmesh;
+   delete V_space;
+   delete W_space;
+   delete M_space;
+   delete AVarf;
+   delete fform;
+   delete dg_coll;
+   delete face;
+
+   if (verbose)
+   {
+      std::cout << "\n\nDone." << std::endl ;
+   }
+
+#ifdef MFEM_USE_PETSC
+   MFEMFinalizePetsc();
+#endif
+   MPI_Finalize();
+
+   return 0;
+}
+
+real_t uFun_ex(const Vector & x)
+{
+   real_t xi(x(0));
+   real_t yi(x(1));
+   int dim = x.Size();
+
+   switch (dim)
+   {
+      case 2:
+      {
+         return 1.0 + xi + sin(2.0*M_PI*xi)*sin(2.0*M_PI*yi);
+         break;
+      }
+      case 3:
+      {
+         real_t zi(x(2));
+         return 1.0*xi + sin(2.0*M_PI*xi)*sin(2.0*M_PI*yi)*sin(2.0*M_PI*zi);
+         break;
+      }
+   }
+
+   return 0;
+}
+
+void qFun_ex(const Vector & x, Vector & q)
+{
+   real_t xi(x(0));
+   real_t yi(x(1));
+   int dim = x.Size();
+
+   switch (dim)
+   {
+      case 2:
+      {
+         q(0) = -diff*1.0 - diff*2.0*M_PI*cos(2.0*M_PI*xi)*sin(2.0*M_PI*yi);
+         q(1) = 0.0 - diff*2.0*M_PI*sin(2.0*M_PI*xi)*cos(2.0*M_PI*yi);
+         break;
+      }
+      case 3:
+      {
+         real_t zi(x(2));
+         q(0) = -diff*1.0 - diff*2.0*M_PI*cos(2.0*M_PI*xi)*sin(2.0*M_PI*yi)*sin(
+                   2.0*M_PI*zi);
+         q(1) = - diff*2.0*M_PI*sin(2.0*M_PI*xi)*cos(2.0*M_PI*yi)*sin(2.0*M_PI*zi);
+         q(2) = - diff*2.0*M_PI*sin(2.0*M_PI*xi)*sin(2.0*M_PI*yi)*cos(2.0*M_PI*zi);
+         break;
+      }
+   }
+}
+
+
+real_t fFun(const Vector & x)
+{
+   real_t xi(x(0));
+   real_t yi(x(1));
+   int dim = x.Size();
+
+   switch (dim)
+   {
+      case 2:
+      {
+         return diff*8.0*M_PI*M_PI*sin(2.0*M_PI*xi)*sin(2.0*M_PI*yi);
+         break;
+      }
+      case 3:
+      {
+         real_t zi(x(2));
+         return diff*12.0*M_PI*M_PI*sin(2.0*M_PI*xi)*sin(2.0*M_PI*yi)*sin(2.0*M_PI*zi);
+         break;
+      }
+   }
+
+   return 0;
+
+}
+
+void pHDGPostProcessing::Postprocessing(ParGridFunction &u_postprocessed)
+{
+   Mesh *mesh = pfes->GetMesh();
+   Array<int>  vdofs;
+   Vector      elmat2, shape, RHS, to_RHS, vals, uvals;
+   real_t      RHS2;
+   DenseMatrix elmat, invdfdx, dshape, dshapedxt, qvals;
+
+   int  ndofs;
+   const FiniteElement *fe_elem;
+   ElementTransformation *Trans;
+
+   for (int i = 0; i < pfes->GetNE(); i++)
+   {
+      pfes->GetElementVDofs(i, vdofs);
+      ndofs = vdofs.Size();
+      vals.SetSize(ndofs);
+      vals = 0.0;
+      elmat.SetSize(ndofs);
+      elmat2.SetSize(ndofs);
+      shape.SetSize(ndofs);
+
+      RHS.SetSize(ndofs);
+      to_RHS.SetSize(ndofs);
+
+      elmat = 0.0;
+      elmat2 = 0.0;
+      RHS = 0.0;
+      RHS2 = 0.0;
+
+      fe_elem = pfes->GetFE(i);
+      int dim = fe_elem->GetDim();
+      int spaceDim = dim;
+      Vector qval_col;
+      qval_col.SetSize(dim);
+      invdfdx.SetSize(dim, spaceDim);
+      dshape.SetSize(ndofs, spaceDim);
+      dshapedxt.SetSize(ndofs, spaceDim);
+
+      Trans = mesh->GetElementTransformation(i);
+
+      const IntegrationRule *ir = IntRule;
+      if (ir == NULL)
+      {
+         int order = 2*fe_elem->GetOrder() + 2;
+         ir = &IntRules.Get(fe_elem->GetGeomType(), order);
+      }
+
+      for (int j = 0; j < ir->GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(j);
+
+         fe_elem->CalcDShape(ip, dshape);
+         fe_elem->CalcShape(ip, shape);
+
+         Trans->SetIntPoint(&ip);
+         // Compute invdfdx = / adj(J),      if J is square
+         //            \ adj(J^t.J).J^t, otherwise
+         CalcAdjugate(Trans->Jacobian(), invdfdx);
+         real_t w = Trans->Weight();
+         w = ip.weight / w;
+         w *= diffcoeff->Eval(*Trans, ip);
+         Mult(dshape, invdfdx, dshapedxt);
+
+         AddMult_a_AAt(w, dshapedxt, elmat);
+
+         dshapedxt *= ip.weight ;
+
+         qval_col = 0.0;
+         for (int ii = 0; ii<dim; ii++)
+         {
+            qval_col(ii) = q->GetValue(i, ip, (ii+1));
+         }
+
+         dshapedxt.Mult(qval_col, to_RHS);
+
+         RHS -= to_RHS;
+
+         shape *= (Trans->Weight() * ip.weight);
+         elmat2 += shape;
+
+         real_t uvalsj;
+         uvalsj = u->GetValue(i, ip, 1);
+
+         real_t rhs_weight = (Trans->Weight() * ip.weight);
+         RHS2  += (uvalsj*rhs_weight);
+
+      }
+
+      // changing the last row and the last entry
+      for (int j = 0; j < ndofs; j++)
+      {
+         elmat(ndofs-1,j) = elmat2(j);
+      }
+      RHS(ndofs-1) = RHS2;
+
+      elmat.Invert();
+      elmat.Mult(RHS, vals);
+      u_postprocessed.SetSubVector(vdofs, vals);
+
+   }
+}

miniapps/hdg/makefile

@@ -0,0 +1,73 @@
+# 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.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../..
+MFEM_BUILD_DIR ?= ../..
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/hdg/,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+# Use the MFEM install directory
+# MFEM_INSTALL_DIR = ../../mfem
+# CONFIG_MK = $(MFEM_INSTALL_DIR)/share/mfem/config.mk
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+SEQ_MINIAPPS = hdg_advection hdg_poisson
+PAR_MINIAPPS = hdg_advectionp hdg_poissonp
+
+ifeq ($(MFEM_USE_MPI),NO)
+   MINIAPPS = $(SEQ_MINIAPPS)
+else
+   MINIAPPS = $(PAR_MINIAPPS) $(SEQ_MINIAPPS)
+endif
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all clean clean-build clean-exec
+
+all: $(MINIAPPS)
+
+# Remove built-in rules
+%: %.cpp
+%.o: %.cpp
+
+# Replace the default implicit rule for *.cpp files
+%: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK) HPP-CPP
+	$(MFEM_CXX) $(MFEM_FLAGS) $< *.o -o $@ $(MFEM_LIBS)
+
+HPP-CPP:
+	$(MFEM_CXX) $(MFEM_FLAGS) -c HDGBilinearForm.cpp
+	$(MFEM_CXX) $(MFEM_FLAGS) -c hdg_integrators.cpp
+
+MFEM_TESTS = MINIAPPS
+include $(MFEM_TEST_MK)
+
+# Testing: Parallel vs. serial runs
+RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
+%-test-par: %
+	@$(call mfem-test,$<, $(RUN_MPI), HDG miniapp)
+%-test-seq: %
+	@$(call mfem-test,$<,, HDG miniapp)
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean: clean-build clean-exec
+
+clean-build:
+	rm -f *.o *~ $(MINIAPPS)
+	rm -rf *.dSYM *.TVD.*breakpoints
+	rm -rf mesh* sol*
+
+clean-exec:
+	@rm -rf *.mesh *.gf