modifications

Mesh reader support for VTK high order elements [feature/vtk_high_order]

CHANGELOG

@@ -35,6 +35,12 @@ Version 4.2.1 (development)
 - Implemented a filter method for the Navier miniapp to stabilize highly
   turbulent flows in direct numerical simulation.
 
+- Added support for reading high-order Lagrange meshes in VTK format. Arbitrary-
+  orders and all element types are supported. See the VTK blog for more info:
+  https://blog.kitware.com/wp-content/uploads/2018/09/Source_Issue_43.pdf
+
+- Added support for reading VTK meshes in XML format.
+
 - Added partial assembly and device support to Example 25/25p, with diagonal
   preconditioning.
 
@@ -52,6 +58,11 @@ Version 4.2.1 (development)
   deprecated EvalSymmetric in MatrixCoefficient. Added DiagonalMatrixCoefficient
   for clarity, which is a typedef of VectorCoefficient.
 
+- Added support for AMG preconditioners for non-symmetric systems (e.g.
+  advection-dominated problems) using hypre's approximate ideal restriction
+  (AIR) AMG. Requires hypre version 2.14.0 or newer. Usage is illustrated in
+  example 9/9p.
+
 - Implemented an adaptive linear solver tolerance option for NewtonSolver based
   on the algorithm of Eisenstat and Walker.
 

NOTICE

@@ -29,4 +29,5 @@ license files. These software products and their licenses are as follows:
 * Catch++ (tests/unit/catch.hpp) -- Boost 1.0 license
 * Gecko (general/gecko.{cpp,hpp}) -- BSD 3-clause license
 * Picojson (fem/picojson.h) -- Custom 2-clause license
+* TinyXML2 (general/tinyxml2.{cpp,h}) -- zlib license
 * Zstr (general/zstr.hpp) -- MIT license

doc/CodeDocumentation.conf.in

@@ -819,7 +819,9 @@ RECURSIVE              = NO
 # run.
 
 EXCLUDE                = @MFEM_SOURCE_DIR@/config/_config.hpp \
-                         @MFEM_SOURCE_DIR@/config/get_hypre_version.cpp
+                         @MFEM_SOURCE_DIR@/config/get_hypre_version.cpp \
+                         @MFEM_SOURCE_DIR@/general/tinyxml2.h \
+                         @MFEM_SOURCE_DIR@/general/tinyxml2.cpp
 
 # The EXCLUDE_SYMLINKS tag can be used to select whether or not files or
 # directories that are symbolic links (a Unix file system feature) are excluded

examples/CMakeLists.txt

@@ -34,6 +34,7 @@ list(APPEND ALL_EXE_SRCS
   ex25.cpp
   ex26.cpp
   ex27.cpp
+  ex91.cpp
   )
 
 if (MFEM_USE_MPI)

examples/ex9.cpp

@@ -290,15 +290,16 @@ int main(int argc, char *argv[])
       k.SetAssemblyLevel(AssemblyLevel::FULL);
    }
    m.AddDomainIntegrator(new MassIntegrator);
-   k.AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
+   constexpr double alpha = -1.0;
+   k.AddDomainIntegrator(new ConvectionIntegrator(velocity, alpha));
    k.AddInteriorFaceIntegrator(
-      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
+      new NonconservativeDGTraceIntegrator(velocity, alpha));
    k.AddBdrFaceIntegrator(
-      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
+      new NonconservativeDGTraceIntegrator(velocity, alpha));
 
    LinearForm b(&fes);
    b.AddBdrFaceIntegrator(
-      new BoundaryFlowIntegrator(inflow, velocity, -1.0, -0.5));
+      new BoundaryFlowIntegrator(inflow, velocity, alpha));
 
    m.Assemble();
    int skip_zeros = 0;

examples/ex91.cpp

@@ -0,0 +1,672 @@
+#include"mfem.hpp"
+
+#include<memory>
+#include<iostream>
+#include<fstream>
+
+namespace mfem {
+
+
+class LinDiffQFunc
+{
+public:
+    LinDiffQFunc(mfem::Coefficient& dd, mfem::Coefficient& ll, double gg_,
+                 double pp0_, double pp1_):diff(dd),load(ll),gg(gg_),pp0(pp0_),pp1(pp1_)
+    {
+
+    }
+
+    double QEnergy(ElementTransformation &T,
+                  const IntegrationPoint &ip,
+                  mfem::Vector& param, mfem::Vector& uu)
+    {
+        double dd=diff.Eval(T,ip);
+        double ll=load.Eval(T,ip);
+
+        double rho0=param[0];
+        double rho1=param[1];
+        double fd=dd*std::pow(rho0,pp0)*std::pow(rho1,pp1);
+
+        double rez = 0.5*(uu[0]*uu[0]+uu[1]*uu[1]+uu[2]*uu[2])*fd
+                        + 0.5*gg*uu[3]*uu[3] -uu[3]*ll;
+
+        return rez;
+    }
+
+    void QResidual(ElementTransformation &T,
+                   const IntegrationPoint &ip,
+                   mfem::Vector& param, mfem::Vector& uu, mfem::Vector& rr)
+    {
+        rr.SetSize(4);
+        double dd=diff.Eval(T,ip);
+        double ll=load.Eval(T,ip);
+
+        double rho0=param[0];
+        double rho1=param[1];
+
+        double fd=dd*std::pow(rho0,pp0)*std::pow(rho1,pp1);
+
+        rr[0]=uu[0]*fd;
+        rr[1]=uu[1]*fd;
+        rr[2]=uu[2]*fd;
+        rr[3]=gg*uu[3]-ll;
+    }
+
+    void AQResidual(ElementTransformation &T,
+                    const IntegrationPoint &ip,
+                    mfem::Vector& param,
+                    mfem::Vector& uu, mfem::Vector& aa, mfem::Vector& rr)
+    {
+         rr.SetSize(2);
+         double dd=diff.Eval(T,ip);
+         double ll=load.Eval(T,ip);
+
+         double rho0=param[0];
+         double rho1=param[1];
+
+         double fd0=dd*pp0*std::pow(rho0,pp0-1.0)*std::pow(rho1,pp1);
+         double fd1=dd*std::pow(rho0,pp0)*pp1*std::pow(rho1,pp1-1.0);
+
+         rr[0] = (aa[0]*uu[0]+aa[1]*uu[1]+aa[2]*uu[2])*fd0;
+         rr[1] = (aa[0]*uu[0]+aa[1]*uu[1]+aa[2]*uu[2])*fd1;
+
+    }
+
+    void QGradResidual(ElementTransformation &T,
+                       const IntegrationPoint &ip,
+                       mfem::Vector& param, mfem::Vector& uu, mfem::DenseMatrix& hh)
+    {
+        hh.SetSize(4);
+        double dd=diff.Eval(T,ip);
+        //double ll=load.Eval(T,ip);
+
+
+        double rho0=param[0];
+        double rho1=param[1];
+
+        double fd=dd*std::pow(rho0,pp0)*std::pow(rho1,pp1);
+        hh=0.0;
+
+        hh(0,0)=fd;
+        hh(1,1)=fd;
+        hh(2,2)=fd;
+        hh(3,3)=gg;
+    }
+
+
+private:
+    mfem::Coefficient& diff;
+    mfem::Coefficient& load;
+    double gg;
+    double pp0;
+    double pp1;
+};
+
+
+class PrmBlockLSFEMDiffusion: public PrmBlockNonlinearFormIntegrator
+{
+public:
+    PrmBlockLSFEMDiffusion(LinDiffQFunc& qfun_)
+    {
+        qfunc=&qfun_;
+    }
+
+    /// Compute the local energy
+    virtual double GetElementEnergy(const Array<const FiniteElement *>&el,
+                                    const Array<const FiniteElement *>&pel,
+                                    ElementTransformation &Tr,
+                                    const Array<const Vector *>&elfun,
+                                    const Array<const Vector *>&pelfun)
+    {
+        int dof_u0 = el[0]->GetDof();
+        int dof_r0 = pel[0]->GetDof();
+        int dof_r1 = pel[1]->GetDof();
+
+        int dim = el[0]->GetDim();
+        int spaceDim = Tr.GetSpaceDim();
+        if (dim != spaceDim)
+        {
+            mfem::mfem_error(" PrmBlockLSFEMDiffusion::GetElementEnergy"
+                             " is not defined on manifold meshes");
+        }
+
+        //shape functions
+        Vector shu0(dof_u0);
+        Vector shr0(dof_r0);
+        Vector shr1(dof_r1);
+        DenseMatrix dsu0(dof_u0,dim);
+        DenseMatrix B(dof_u0, 4);
+        B=0.0;
+
+        double w;
+
+        Vector param(2); param=0.0;
+        Vector uu(4); uu=0.0;
+
+        double energy =0.0;
+
+        const IntegrationRule *ir = nullptr;
+        if(ir==nullptr){
+            int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
+                                +pel[0]->GetOrder()+pel[1]->GetOrder();
+            ir=&IntRules.Get(Tr.GetGeometryType(),order);
+        }
+
+        for (int i = 0; i < ir->GetNPoints(); i++)
+        {
+            const IntegrationPoint &ip = ir->IntPoint(i);
+            Tr.SetIntPoint(&ip);
+            w=Tr.Weight();
+            w = ip.weight * w;
+
+            el[0]->CalcPhysDShape(Tr,dsu0);
+            el[0]->CalcPhysShape(Tr,shu0);
+            pel[0]->CalcPhysShape(Tr,shr0);
+            pel[1]->CalcPhysShape(Tr,shr1);
+
+            param[0]=shr0*(*pelfun[0]);
+            param[1]=shr1*(*pelfun[1]);
+
+            //set the matrix B
+            for(int jj=0;jj<dim;jj++)
+            {
+                B.SetCol(jj,dsu0.GetColumn(jj));
+            }
+            B.SetCol(3,shu0);
+            B.MultTranspose(*elfun[0],uu);
+            energy=energy+w * qfunc->QEnergy(Tr,ip,param,uu);
+        }
+        return energy;
+    }
+
+    /// Perform the local action of the BlockNonlinearFormIntegrator
+    virtual void AssembleElementVector(const Array<const FiniteElement *> &el,
+                                       const Array<const FiniteElement *>&pel,
+                                       ElementTransformation &Tr,
+                                       const Array<const Vector *> &elfun,
+                                       const Array<const Vector *>&pelfun,
+                                       const Array<Vector *> &elvec)
+    {
+        int dof_u0 = el[0]->GetDof();
+        int dof_r0 = pel[0]->GetDof();
+        int dof_r1 = pel[1]->GetDof();
+
+        int dim = el[0]->GetDim();
+
+        elvec[0]->SetSize(dof_u0);
+        *elvec[0]=0.0;
+        int spaceDim = Tr.GetSpaceDim();
+        if (dim != spaceDim)
+        {
+            mfem::mfem_error(" PrmBlockLSFEMDiffusion::AssembleElementVector"
+                             " is not defined on manifold meshes");
+        }
+
+        //shape functions
+        Vector shu0(dof_u0);
+        Vector shr0(dof_r0);
+        Vector shr1(dof_r1);
+        DenseMatrix dsu0(dof_u0,dim);
+        DenseMatrix B(dof_u0, 4);
+        B=0.0;
+
+        double w;
+
+        Vector param(2); param=0.0;
+        Vector uu(4); uu=0.0;
+        Vector rr;
+        Vector lvec; lvec.SetSize(dof_u0);
+
+        const IntegrationRule *ir = nullptr;
+        if(ir==nullptr){
+            int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
+                                +pel[0]->GetOrder()+pel[1]->GetOrder();
+            ir=&IntRules.Get(Tr.GetGeometryType(),order);
+        }
+
+        for (int i = 0; i < ir->GetNPoints(); i++)
+        {
+            const IntegrationPoint &ip = ir->IntPoint(i);
+            Tr.SetIntPoint(&ip);
+            w=Tr.Weight();
+            w = ip.weight * w;
+
+            el[0]->CalcPhysDShape(Tr,dsu0);
+            el[0]->CalcPhysShape(Tr,shu0);
+            pel[0]->CalcPhysShape(Tr,shr0);
+            pel[1]->CalcPhysShape(Tr,shr1);
+
+            param[0]=shr0*(*pelfun[0]);
+            param[1]=shr1*(*pelfun[1]);
+
+            //set the matrix B
+            for(int jj=0;jj<dim;jj++)
+            {
+                B.SetCol(jj,dsu0.GetColumn(jj));
+            }
+            B.SetCol(3,shu0);
+            B.MultTranspose(*elfun[0],uu);
+            qfunc->QResidual(Tr,ip,param, uu, rr);
+
+            B.Mult(rr,lvec);
+            elvec[0]->Add(w,lvec);
+        }
+
+    }
+
+    virtual void AssembleFaceVector(const Array<const FiniteElement *> &el1,
+                                    const Array<const FiniteElement *> &el2,
+                                    const Array<const FiniteElement *> &pel1,
+                                    const Array<const FiniteElement *> &pel2,
+                                    FaceElementTransformations &Tr,
+                                    const Array<const Vector *> &elfun,
+                                    const Array<const Vector *>&pelfun,
+                                    const Array<Vector *> &elvect)
+    {
+
+    }
+
+    /// Assemble the local gradient matrix
+    virtual void AssembleElementGrad(const Array<const FiniteElement*> &el,
+                                     const Array<const FiniteElement *>&pel,
+                                     ElementTransformation &Tr,
+                                     const Array<const Vector *> &elfun,
+                                     const Array<const Vector *>&pelfun,
+                                     const Array2D<DenseMatrix *> &elmats)
+    {
+        int dof_u0 = el[0]->GetDof();
+        int dof_r0 = pel[0]->GetDof();
+        int dof_r1 = pel[1]->GetDof();
+
+        int dim = el[0]->GetDim();
+
+        //elmats[0]->Size(dof_u0, dof_u0);
+        //*elmats[0]=0.0;
+
+        DenseMatrix* K=elmats(0,0);
+        K->SetSize(dof_u0,dof_u0);
+        (*K)=0.0;
+
+        int spaceDim = Tr.GetSpaceDim();
+        if (dim != spaceDim)
+        {
+            mfem::mfem_error(" PrmBlockLSFEMDiffusion::AssembleElementVector"
+                             " is not defined on manifold meshes");
+        }
+
+        //shape functions
+        Vector shu0(dof_u0);
+        Vector shr0(dof_r0);
+        Vector shr1(dof_r1);
+        DenseMatrix dsu0(dof_u0,dim);
+        DenseMatrix B(dof_u0, 4);
+        DenseMatrix A(dof_u0, 4);
+        B=0.0;
+
+        double w;
+
+        Vector param(2); param=0.0;
+        Vector uu(4); uu=0.0;
+        DenseMatrix hh;
+        Vector lvec; lvec.SetSize(dof_u0);
+
+        const IntegrationRule *ir = nullptr;
+        if(ir==nullptr){
+            int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
+                                +pel[0]->GetOrder()+pel[1]->GetOrder();
+            ir=&IntRules.Get(Tr.GetGeometryType(),order);
+        }
+
+        for (int i = 0; i < ir->GetNPoints(); i++)
+        {
+            const IntegrationPoint &ip = ir->IntPoint(i);
+            Tr.SetIntPoint(&ip);
+            w = Tr.Weight();
+            w = ip.weight * w;
+
+            el[0]->CalcPhysDShape(Tr,dsu0);
+            el[0]->CalcPhysShape(Tr,shu0);
+            pel[0]->CalcPhysShape(Tr,shr0);
+            pel[1]->CalcPhysShape(Tr,shr1);
+
+            param[0]=shr0*(*pelfun[0]);
+            param[1]=shr1*(*pelfun[1]);
+
+            //set the matrix B
+            for(int jj=0;jj<dim;jj++)
+            {
+                B.SetCol(jj,dsu0.GetColumn(jj));
+            }
+            B.SetCol(3,shu0);
+            B.MultTranspose(*elfun[0],uu);
+            qfunc->QGradResidual(Tr,ip,param,uu,hh);
+            Mult(B,hh,A);
+            AddMult_a_ABt(w,A,B,*K);
+        }
+    }
+
+    virtual void AssembleFaceGrad(const Array<const FiniteElement *>&el1,
+                                  const Array<const FiniteElement *>&el2,
+                                  const Array<const FiniteElement *> &pel1,
+                                  const Array<const FiniteElement *> &pel2,
+                                  FaceElementTransformations &Tr,
+                                  const Array<const Vector *> &elfun,
+                                  const Array<const Vector *>&pelfun,
+                                  const Array2D<DenseMatrix *> &elmats)
+    {
+
+    }
+
+    virtual void AssemblePrmElementVector(const Array<const FiniteElement *> &el,
+                                          const Array<const FiniteElement *> &pel,
+                                          ElementTransformation &Tr,
+                                          const Array<const Vector *> &elfun,
+                                          const Array<const Vector *> &alfun,
+                                          const Array<const Vector *> &pelfun,
+                                          const Array<Vector *> &elvec)
+    {
+        int dof_u0 = el[0]->GetDof();
+        int dof_r0 = pel[0]->GetDof();
+        int dof_r1 = pel[1]->GetDof();
+
+        int dim = el[0]->GetDim();
+
+        Vector& e0 = *(elvec[0]);
+        Vector& e1 = *(elvec[1]);
+
+        e0.SetSize(dof_r0);
+        e0=0.0;
+        e1.SetSize(dof_r1);
+        e1=0.0;
+
+        int spaceDim = Tr.GetSpaceDim();
+        if (dim != spaceDim)
+        {
+            mfem::mfem_error(" PrmBlockLSFEMDiffusion::AssembleElementVector"
+                             " is not defined on manifold meshes");
+        }
+
+        //shape functions
+        Vector shu0(dof_u0);
+        Vector shr0(dof_r0);
+        Vector shr1(dof_r1);
+        DenseMatrix dsu0(dof_u0,dim);
+        DenseMatrix B(dof_u0, 4);
+        B=0.0;
+
+        double w;
+
+        Vector param(2); param=0.0;
+        Vector uu(4); uu=0.0;
+        Vector aa(4); aa=0.0;
+        Vector rr;
+        Vector lvec0; lvec0.SetSize(dof_r0);
+        Vector lvec1; lvec1.SetSize(dof_r1);
+
+        const IntegrationRule *ir = nullptr;
+        if(ir==nullptr){
+            int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
+                                +pel[0]->GetOrder()+pel[1]->GetOrder();
+            ir=&IntRules.Get(Tr.GetGeometryType(),order);
+        }
+
+        for (int i = 0; i < ir->GetNPoints(); i++)
+        {
+            const IntegrationPoint &ip = ir->IntPoint(i);
+            Tr.SetIntPoint(&ip);
+            w=Tr.Weight();
+            w = ip.weight * w;
+
+            el[0]->CalcPhysDShape(Tr,dsu0);
+            el[0]->CalcPhysShape(Tr,shu0);
+            pel[0]->CalcPhysShape(Tr,shr0);
+            pel[1]->CalcPhysShape(Tr,shr1);
+
+            param[0]=shr0*(*pelfun[0]);
+            param[1]=shr1*(*pelfun[1]);
+
+            //set the matrix B
+            for(int jj=0;jj<dim;jj++)
+            {
+                B.SetCol(jj,dsu0.GetColumn(jj));
+            }
+            B.SetCol(3,shu0);
+            B.MultTranspose(*elfun[0],uu);
+            B.MultTranspose(*alfun[0],aa);
+
+            qfunc->AQResidual(Tr, ip, param, uu, aa, rr);
+
+            lvec0=shr0;
+            lvec0*=rr[0];
+            lvec1=shr1;
+            lvec1*=rr[1];
+
+            e0.Add(w,lvec0);
+            e1.Add(w,lvec1);
+        }
+
+    }
+
+    virtual void AssemblePrmFaceVector(const Array<const FiniteElement *> &el1,
+                                       const Array<const FiniteElement *> &el2,
+                                       const Array<const FiniteElement *> &pel1,
+                                       const Array<const FiniteElement *> &pel2,
+                                       FaceElementTransformations &Tr,
+                                       const Array<const Vector *> &elfun,
+                                       const Array<const Vector *> &alfun,
+                                       const Array<const Vector *> &pelfun,
+                                       const Array<Vector *> &elvect)
+    {
+
+    }
+
+
+
+private:
+    LinDiffQFunc* qfunc;
+};
+
+}
+
+int main(int argc, char *argv[])
+{
+    const char *mesh_file = "../../data/beam-tet.mesh";
+    int ser_ref_levels = 1;
+    int order = 2;
+    bool visualization = true;
+    double newton_rel_tol = 1e-4;
+    double newton_abs_tol = 1e-6;
+    int newton_iter = 10;
+    int print_level = 0;
+
+    mfem::OptionsParser args(argc, argv);
+    args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+    args.AddOption(&ser_ref_levels,
+                     "-rs",
+                     "--refine-serial",
+                     "Number of times to refine the mesh uniformly in serial.");
+    args.AddOption(&order,
+                      "-o",
+                      "--order",
+                      "Order (degree) of the finite elements.");
+    args.AddOption(&visualization,
+                      "-vis",
+                      "--visualization",
+                      "-no-vis",
+                      "--no-visualization",
+                      "Enable or disable GLVis visualization.");
+    args.AddOption(&newton_rel_tol,
+                      "-rel",
+                      "--relative-tolerance",
+                      "Relative tolerance for the Newton solve.");
+    args.AddOption(&newton_abs_tol,
+                      "-abs",
+                      "--absolute-tolerance",
+                      "Absolute tolerance for the Newton solve.");
+    args.AddOption(&newton_iter,
+                      "-it",
+                      "--newton-iterations",
+                      "Maximum iterations for the Newton solve.");
+
+    args.Parse();
+    if (!args.Good())
+    {
+        args.PrintUsage(std::cout);
+        return 1;
+    }
+    args.PrintOptions(std::cout);
+    // 3. Read the (serial) mesh from the given mesh file on all processors.  We
+    //    can handle triangular, quadrilateral, tetrahedral and hexahedral meshes
+    //    with the same code.
+    mfem::Mesh *mesh = new mfem::Mesh(mesh_file, 1, 1);
+    int dim = mesh->Dimension();
+
+    // 4. Refine the mesh in serial to increase the resolution. In this example
+    //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
+    //    a command-line parameter.
+    for (int lev = 0; lev < ser_ref_levels; lev++)
+    {
+        mesh->UniformRefinement();
+    }
+
+    /// Define the q-function
+    mfem::ConstantCoefficient* dc=new mfem::ConstantCoefficient(1.0);
+    mfem::ConstantCoefficient* lc=new mfem::ConstantCoefficient(1.0);
+    mfem::LinDiffQFunc* qfun=new mfem::LinDiffQFunc(*dc,*lc,1.0,1.0,1.0);
+
+    mfem::H1_FECollection fec00(order, dim);
+    mfem::L2_FECollection fec01(order, dim);
+    mfem::FiniteElementSpace* bfes00=new mfem::FiniteElementSpace(mesh,&fec00,1,mfem::Ordering::byVDIM);
+    mfem::FiniteElementSpace* pfes00=new mfem::FiniteElementSpace(mesh,&fec00,1,mfem::Ordering::byVDIM);
+    mfem::FiniteElementSpace* pfes01=new mfem::FiniteElementSpace(mesh,&fec01,1,mfem::Ordering::byVDIM);
+
+    /// Define parametric nonlinear form
+    mfem::Array<mfem::FiniteElementSpace*> bfes;
+    mfem::Array<mfem::FiniteElementSpace*> pfes;
+
+    bfes.Append(bfes00);
+    pfes.Append(pfes00);
+    pfes.Append(pfes01);
+
+    mfem::PrmBlockNonlinearForm* nf=new mfem::PrmBlockNonlinearForm(bfes,pfes);
+    nf->AddDomainIntegrator(new mfem::PrmBlockLSFEMDiffusion(*qfun));
+
+    /// Define the grid functions
+    mfem::GridFunction* bgf00=new mfem::GridFunction(bfes00);
+    mfem::GridFunction* pgf00=new mfem::GridFunction(pfes00);
+    mfem::GridFunction* pgf01=new mfem::GridFunction(pfes01);
+    mfem::GridFunction* ggf00=new mfem::GridFunction(pfes00);
+    mfem::GridFunction* ggf01=new mfem::GridFunction(pfes01);
+
+    *bgf00=0.0;
+    *pgf00=1.0;
+    *pgf01=1.0;
+
+    mfem::BlockVector solbv; solbv.Update(nf->GetBlockTrueOffsets()); solbv=0.0;
+    mfem::BlockVector resbv; resbv.Update(nf->GetBlockTrueOffsets()); resbv=0.0;
+    mfem::BlockVector adjbv; adjbv.Update(nf->GetBlockTrueOffsets()); adjbv=0.0;
+    mfem::BlockVector prmbv; prmbv.Update(nf->PrmGetBlockTrueOffsets()); prmbv=1.0;
+    mfem::BlockVector grdbv; grdbv.Update(nf->PrmGetBlockTrueOffsets()); grdbv=0.0;
+
+    bgf00->SetFromTrueDofs(solbv.GetBlock(0));
+    pgf00->SetFromTrueDofs(prmbv.GetBlock(0));
+    pgf01->SetFromTrueDofs(prmbv.GetBlock(1));
+
+    nf->SetPrmFields(prmbv);
+    double energy = nf->GetEnergy(solbv);
+
+    nf->Mult(solbv,resbv);
+    std::cout<<"Norm res="<<resbv.Norml2()<<std::endl;
+
+    //mfem::Operator& K=nf->GetGradient(solbv);
+    std::cout<<"energy ="<< energy<<std::endl;
+
+    nf->SetStateFields(solbv);
+    nf->SetAdjointFields(adjbv);
+    nf->PrmMult(prmbv,grdbv);
+
+
+    //set the BC for the physics
+    mfem::Array<mfem::Array<int> *> ess_bdr;
+    mfem::Array<mfem::Vector*>      ess_rhs;
+    ess_bdr.Append(new mfem::Array<int>(mesh->bdr_attributes.Max()));
+    ess_rhs.Append(nullptr);
+    (*ess_bdr[0]) = 1;
+    nf->SetEssentialBC(ess_bdr,ess_rhs);
+
+    //define the solvers
+    mfem::UMFPackSolver* umfsolv=new mfem::UMFPackSolver();
+
+    mfem::GMRESSolver *gmres;
+    gmres = new mfem::GMRESSolver();
+    gmres->SetAbsTol(newton_abs_tol/10);
+    gmres->SetRelTol(newton_rel_tol/10);
+    gmres->SetMaxIter(100);
+    gmres->SetPrintLevel(print_level);
+    //gmres->SetPreconditioner(*prec);
+
+
+    mfem::NewtonSolver *ns;
+    ns = new mfem::NewtonSolver();
+    ns->iterative_mode = true;
+    ns->SetSolver(*gmres);
+    ns->SetOperator(*nf);
+    ns->SetPrintLevel(print_level);
+    ns->SetRelTol(newton_rel_tol);
+    ns->SetAbsTol(newton_abs_tol);
+    ns->SetMaxIter(newton_iter);
+
+    mfem::Vector b; //RHS is zero
+    solbv=0.0;
+    ns->Mult(b, solbv);
+
+
+    nf->SetStateFields(solbv);
+    nf->SetAdjointFields(solbv);
+    nf->PrmMult(prmbv,grdbv);
+
+
+
+
+    mfem::ParaViewDataCollection *dacol = new mfem::ParaViewDataCollection("Example91",
+                                                                      mesh);
+
+    ggf00->SetFromTrueDofs(grdbv.GetBlock(0));
+    ggf01->SetFromTrueDofs(grdbv.GetBlock(1));
+    pgf00->SetFromTrueDofs(solbv.GetBlock(0));
+
+    dacol->SetLevelsOfDetail(order);
+    dacol->RegisterField("sol", pgf00);
+    dacol->RegisterField("grad00", ggf00);
+    dacol->RegisterField("grad01", ggf01);
+
+    dacol->SetTime(1.0);
+    dacol->SetCycle(1);
+    dacol->Save();
+
+    delete dacol;
+
+    delete ns;
+    delete umfsolv;
+    delete gmres;
+    delete ess_bdr[0];
+
+    delete bgf00;
+    delete pgf00;
+    delete pgf01;
+    delete ggf00;
+    delete ggf01;
+
+
+    delete nf;
+    delete pfes01;
+    delete pfes00;
+    delete bfes00;
+
+    delete qfun;
+    delete lc;
+    delete dc;
+
+    delete mesh;
+
+}

examples/ex9p.cpp

@@ -64,6 +64,66 @@ double inflow_function(const Vector &x);
 // Mesh bounding box
 Vector bb_min, bb_max;
 
+// Type of preconditioner for implicit time integrator
+enum class PrecType : int
+{
+   ILU = 0,
+   AIR = 1
+};
+
+#if MFEM_HYPRE_VERSION >= 21800
+// Algebraic multigrid preconditioner for advective problems based on
+// approximate ideal restriction (AIR). Most effective when matrix is
+// first scaled by DG block inverse, and AIR applied to scaled matrix.
+// See https://doi.org/10.1137/17M1144350.
+class AIR_prec : public Solver
+{
+private:
+   const HypreParMatrix *A;
+   // Copy of A scaled by block-diagonal inverse
+   HypreParMatrix A_s;
+
+   HypreBoomerAMG *AIR_solver;
+   int blocksize;
+
+public:
+   AIR_prec(int blocksize_) : AIR_solver(NULL), blocksize(blocksize_) { }
+
+   void SetOperator(const Operator &op)
+   {
+      width = op.Width();
+      height = op.Height();
+
+      A = dynamic_cast<const HypreParMatrix *>(&op);
+      MFEM_VERIFY(A != NULL, "AIR_prec requires a HypreParMatrix.")
+
+      // Scale A by block-diagonal inverse
+      BlockInverseScale(A, &A_s, NULL, NULL, blocksize,
+                        BlockInverseScaleJob::MATRIX_ONLY);
+      delete AIR_solver;
+      AIR_solver = new HypreBoomerAMG(A_s);
+      AIR_solver->SetAdvectiveOptions(1, "", "FA");
+      AIR_solver->SetPrintLevel(0);
+      AIR_solver->SetMaxLevels(50);
+   }
+
+   virtual void Mult(const Vector &x, Vector &y) const
+   {
+      // Scale the rhs by block inverse and solve system
+      HypreParVector z_s;
+      BlockInverseScale(A, NULL, &x, &z_s, blocksize,
+                        BlockInverseScaleJob::RHS_ONLY);
+      AIR_solver->Mult(z_s, y);
+   }
+
+   ~AIR_prec()
+   {
+      delete AIR_solver;
+   }
+};
+#endif
+
+
 class DG_Solver : public Solver
 {
 private:
@@ -71,24 +131,37 @@ private:
    SparseMatrix M_diag;
    HypreParMatrix *A;
    GMRESSolver linear_solver;
-   BlockILU prec;
+   Solver *prec;
    double dt;
 public:
-   DG_Solver(HypreParMatrix &M_, HypreParMatrix &K_, const FiniteElementSpace &fes)
+   DG_Solver(HypreParMatrix &M_, HypreParMatrix &K_, const FiniteElementSpace &fes,
+             PrecType prec_type)
       : M(M_),
         K(K_),
         A(NULL),
         linear_solver(M.GetComm()),
-        prec(fes.GetFE(0)->GetDof(),
-             BlockILU::Reordering::MINIMUM_DISCARDED_FILL),
         dt(-1.0)
    {
+      int block_size = fes.GetFE(0)->GetDof();
+      if (prec_type == PrecType::ILU)
+      {
+         prec = new BlockILU(block_size,
+                             BlockILU::Reordering::MINIMUM_DISCARDED_FILL);
+      }
+      else if (prec_type == PrecType::AIR)
+      {
+#if MFEM_HYPRE_VERSION >= 21800
+         prec = new AIR_prec(block_size);
+#else
+         MFEM_ABORT("Must have MFEM_HYPRE_VERSION >= 21800 to use AIR.\n");
+#endif
+      }
       linear_solver.iterative_mode = false;
       linear_solver.SetRelTol(1e-9);
       linear_solver.SetAbsTol(0.0);
       linear_solver.SetMaxIter(100);
       linear_solver.SetPrintLevel(0);
-      linear_solver.SetPreconditioner(prec);
+      linear_solver.SetPreconditioner(*prec);
 
       M.GetDiag(M_diag);
    }
@@ -121,10 +194,12 @@ public:
 
    ~DG_Solver()
    {
+      delete prec;
       delete A;
    }
 };
 
+
 /** A time-dependent operator for the right-hand side of the ODE. The DG weak
     form of du/dt = -v.grad(u) is M du/dt = K u + b, where M and K are the mass
     and advection matrices, and b describes the flow on the boundary. This can
@@ -142,7 +217,8 @@ private:
    mutable Vector z;
 
 public:
-   FE_Evolution(ParBilinearForm &_M, ParBilinearForm &_K, const Vector &_b);
+   FE_Evolution(ParBilinearForm &_M, ParBilinearForm &_K, const Vector &_b,
+                PrecType prec_type);
 
    virtual void Mult(const Vector &x, Vector &y) const;
    virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);
@@ -154,10 +230,9 @@ public:
 int main(int argc, char *argv[])
 {
    // 1. Initialize MPI.
-   int num_procs, myid;
-   MPI_Init(&argc, &argv);
-   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
-   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+   MPI_Session mpi;
+   int num_procs = mpi.WorldSize();
+   int myid = mpi.WorldRank();
 
    // 2. Parse command-line options.
    problem = 0;
@@ -178,7 +253,11 @@ int main(int argc, char *argv[])
    bool adios2 = false;
    bool binary = false;
    int vis_steps = 5;
-
+#if MFEM_HYPRE_VERSION >= 21800
+   PrecType prec_type = PrecType::AIR;
+#else
+   PrecType prec_type = PrecType::ILU;
+#endif
    int precision = 8;
    cout.precision(precision);
 
@@ -212,6 +291,8 @@ int main(int argc, char *argv[])
                   "Final time; start time is 0.");
    args.AddOption(&dt, "-dt", "--time-step",
                   "Time step.");
+   args.AddOption((int *)&prec_type, "-pt", "--prec-type", "Preconditioner for "
+                  "implicit solves. 0 for ILU, 1 for pAIR-AMG.");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
@@ -232,20 +313,19 @@ int main(int argc, char *argv[])
    args.Parse();
    if (!args.Good())
    {
-      if (myid == 0)
+      if (mpi.Root())
       {
          args.PrintUsage(cout);
       }
-      MPI_Finalize();
       return 1;
    }
-   if (myid == 0)
+   if (mpi.Root())
    {
       args.PrintOptions(cout);
    }
 
    Device device(device_config);
-   if (myid == 0) { device.Print(); }
+   if (mpi.Root()) { device.Print(); }
 
    // 3. Read the serial mesh from the given mesh file on all processors. We can
    //    handle geometrically periodic meshes in this code.
@@ -272,12 +352,11 @@ int main(int argc, char *argv[])
       case 23: ode_solver = new SDIRK23Solver; break;
       case 24: ode_solver = new SDIRK34Solver; break;
       default:
-         if (myid == 0)
+         if (mpi.Root())
          {
             cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
          }
          delete mesh;
-         MPI_Finalize();
          return 3;
    }
 
@@ -311,7 +390,7 @@ int main(int argc, char *argv[])
    ParFiniteElementSpace *fes = new ParFiniteElementSpace(pmesh, &fec);
 
    HYPRE_Int global_vSize = fes->GlobalTrueVSize();
-   if (myid == 0)
+   if (mpi.Root())
    {
       cout << "Number of unknowns: " << global_vSize << endl;
    }
@@ -342,15 +421,16 @@ int main(int argc, char *argv[])
    }
 
    m->AddDomainIntegrator(new MassIntegrator);
-   k->AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
+   constexpr double alpha = -1.0;
+   k->AddDomainIntegrator(new ConvectionIntegrator(velocity, alpha));
    k->AddInteriorFaceIntegrator(
-      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
+      new NonconservativeDGTraceIntegrator(velocity, alpha));
    k->AddBdrFaceIntegrator(
-      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
+      new NonconservativeDGTraceIntegrator(velocity, alpha));
 
    ParLinearForm *b = new ParLinearForm(fes);
    b->AddBdrFaceIntegrator(
-      new BoundaryFlowIntegrator(inflow, velocity, -1.0, -0.5));
+      new BoundaryFlowIntegrator(inflow, velocity, alpha));
 
    int skip_zeros = 0;
    m->Assemble();
@@ -451,11 +531,11 @@ int main(int argc, char *argv[])
       sout.open(vishost, visport);
       if (!sout)
       {
-         if (myid == 0)
+         if (mpi.Root())
             cout << "Unable to connect to GLVis server at "
                  << vishost << ':' << visport << endl;
          visualization = false;
-         if (myid == 0)
+         if (mpi.Root())
          {
             cout << "GLVis visualization disabled.\n";
          }
@@ -467,7 +547,7 @@ int main(int argc, char *argv[])
          sout << "solution\n" << *pmesh << *u;
          sout << "pause\n";
          sout << flush;
-         if (myid == 0)
+         if (mpi.Root())
             cout << "GLVis visualization paused."
                  << " Press space (in the GLVis window) to resume it.\n";
       }
@@ -476,7 +556,7 @@ int main(int argc, char *argv[])
    // 10. Define the time-dependent evolution operator describing the ODE
    //     right-hand side, and perform time-integration (looping over the time
    //     iterations, ti, with a time-step dt).
-   FE_Evolution adv(*m, *k, *B);
+   FE_Evolution adv(*m, *k, *B, prec_type);
 
    double t = 0.0;
    adv.SetTime(t);
@@ -493,7 +573,7 @@ int main(int argc, char *argv[])
 
       if (done || ti % vis_steps == 0)
       {
-         if (myid == 0)
+         if (mpi.Root())
          {
             cout << "time step: " << ti << ", time: " << t << endl;
          }
@@ -564,16 +644,14 @@ int main(int argc, char *argv[])
 #endif
    delete dc;
 
-   MPI_Finalize();
    return 0;
 }
 
 
 // Implementation of class FE_Evolution
 FE_Evolution::FE_Evolution(ParBilinearForm &_M, ParBilinearForm &_K,
-                           const Vector &_b)
-   : TimeDependentOperator(_M.Height()),
-     b(_b),
+                           const Vector &_b, PrecType prec_type)
+   : TimeDependentOperator(_M.Height()), b(_b),
      M_solver(_M.ParFESpace()->GetComm()),
      z(_M.Height())
 {
@@ -598,7 +676,7 @@ FE_Evolution::FE_Evolution(ParBilinearForm &_M, ParBilinearForm &_K,
       HypreSmoother *hypre_prec = new HypreSmoother(M_mat, HypreSmoother::Jacobi);
       M_prec = hypre_prec;
 
-      dg_solver = new DG_Solver(M_mat, K_mat, *_M.FESpace());
+      dg_solver = new DG_Solver(M_mat, K_mat, *_M.FESpace(), prec_type);
    }
    else
    {
@@ -614,6 +692,10 @@ FE_Evolution::FE_Evolution(ParBilinearForm &_M, ParBilinearForm &_K,
    M_solver.SetPrintLevel(0);
 }
 
+// Solve the equation:
+//    u_t = M^{-1}(Ku + b),
+// by solving associated linear system
+//    (M - dt*K) d = K*u + b
 void FE_Evolution::ImplicitSolve(const double dt, const Vector &x, Vector &k)
 {
    K->Mult(x, z);

fem/CMakeLists.txt

@@ -66,6 +66,7 @@ set(SRCS
   tmop_tools.cpp
   gslib.cpp
   transfer.cpp
+  prmnonlinearform.cpp
   )
 
 set(HDRS
@@ -112,6 +113,7 @@ set(HDRS
   tmop_tools.hpp
   gslib.hpp
   transfer.hpp
+  prmnonlinearform.hpp
   )
 
 if (MFEM_USE_SIDRE)
@@ -136,7 +138,8 @@ if (MFEM_USE_MPI)
     pgridfunc.cpp
     plinearform.cpp
     pnonlinearform.cpp
-    prestriction.cpp)
+    prestriction.cpp
+    pprmnonlinearform.cpp)
   # If this list (HDRS -> HEADERS) is used for install, we probably want the
   # headers added all the time.
   list(APPEND HDRS
@@ -145,7 +148,8 @@ if (MFEM_USE_MPI)
     pgridfunc.hpp
     plinearform.hpp
     pnonlinearform.hpp
-    prestriction.hpp)
+    prestriction.hpp
+    pprmnonlinearform.hpp)
 endif()
 
 convert_filenames_to_full_paths(SRCS)

fem/bilinearform.hpp

@@ -854,9 +854,9 @@ public:
 
       This returns the same operator as FormRectangularLinearSystem(), but does
       without the transformations of the right-hand side. */
-   void FormRectangularSystemMatrix(const Array<int> &trial_tdof_list,
-                                    const Array<int> &test_tdof_list,
-                                    OperatorHandle &A);
+   virtual void FormRectangularSystemMatrix(const Array<int> &trial_tdof_list,
+                                            const Array<int> &test_tdof_list,
+                                            OperatorHandle &A);
 
    /** @brief Form the column-constrained linear system matrix A.
        See FormRectangularSystemMatrix() for details.
@@ -883,10 +883,11 @@ public:
        Return in @a A a *reference* to the system matrix that is column-constrained.
        The reference will be invalidated when SetOperatorType(), Update(), or the
        destructor is called. */
-   void FormRectangularLinearSystem(const Array<int> &trial_tdof_list,
-                                    const Array<int> &test_tdof_list,
-                                    Vector &x, Vector &b,
-                                    OperatorHandle &A, Vector &X, Vector &B);
+   virtual void FormRectangularLinearSystem(const Array<int> &trial_tdof_list,
+                                            const Array<int> &test_tdof_list,
+                                            Vector &x, Vector &b,
+                                            OperatorHandle &A, Vector &X,
+                                            Vector &B);
 
    /** @brief Form the linear system A X = B, corresponding to this bilinear
        form and the linear form @a b(.).

fem/bilininteg.cpp

@@ -229,6 +229,159 @@ void SumIntegrator::AssembleElementMatrix(
    }
 }
 
+void SumIntegrator::AssembleElementMatrix2(
+   const FiniteElement &el1, const FiniteElement &el2,
+   ElementTransformation &Trans, DenseMatrix &elmat)
+{
+   MFEM_ASSERT(integrators.Size() > 0, "empty SumIntegrator.");
+
+   integrators[0]->AssembleElementMatrix2(el1, el2, Trans, elmat);
+   for (int i = 1; i < integrators.Size(); i++)
+   {
+      integrators[i]->AssembleElementMatrix2(el1, el2, Trans, elem_mat);
+      elmat += elem_mat;
+   }
+}
+
+void SumIntegrator::AssembleFaceMatrix(
+   const FiniteElement &el1, const FiniteElement &el2,
+   FaceElementTransformations &Trans, DenseMatrix &elmat)
+{
+   MFEM_ASSERT(integrators.Size() > 0, "empty SumIntegrator.");
+
+   integrators[0]->AssembleFaceMatrix(el1, el2, Trans, elmat);
+   for (int i = 1; i < integrators.Size(); i++)
+   {
+      integrators[i]->AssembleFaceMatrix(el1, el2, Trans, elem_mat);
+      elmat += elem_mat;
+   }
+}
+
+void SumIntegrator::AssembleFaceMatrix(
+   const FiniteElement &tr_fe,
+   const FiniteElement &te_fe1, const FiniteElement &te_fe2,
+   FaceElementTransformations &Trans, DenseMatrix &elmat)
+{
+   MFEM_ASSERT(integrators.Size() > 0, "empty SumIntegrator.");
+
+   integrators[0]->AssembleFaceMatrix(tr_fe, te_fe1, te_fe2, Trans, elmat);
+   for (int i = 1; i < integrators.Size(); i++)
+   {
+      integrators[i]->AssembleFaceMatrix(tr_fe, te_fe1, te_fe2, Trans, elem_mat);
+      elmat += elem_mat;
+   }
+}
+
+void SumIntegrator::AssemblePA(const FiniteElementSpace& fes)
+{
+   for (int i = 0; i < integrators.Size(); i++)
+   {
+      integrators[i]->AssemblePA(fes);
+   }
+}
+
+void SumIntegrator::AssembleDiagonalPA(Vector &diag)
+{
+   for (int i = 0; i < integrators.Size(); i++)
+   {
+      integrators[i]->AssembleDiagonalPA(diag);
+   }
+}
+
+void SumIntegrator::AssemblePAInteriorFaces(const FiniteElementSpace &fes)
+{
+   for (int i = 0; i < integrators.Size(); i++)
+   {
+      integrators[i]->AssemblePAInteriorFaces(fes);
+   }
+}
+
+void SumIntegrator::AssemblePABoundaryFaces(const FiniteElementSpace &fes)
+{
+   for (int i = 0; i < integrators.Size(); i++)
+   {
+      integrators[i]->AssemblePABoundaryFaces(fes);
+   }
+}
+
+void SumIntegrator::AddMultPA(const Vector& x, Vector& y) const
+{
+   for (int i = 0; i < integrators.Size(); i++)
+   {
+      integrators[i]->AddMultPA(x, y);
+   }
+}
+
+void SumIntegrator::AddMultTransposePA(const Vector &x, Vector &y) const
+{
+   for (int i = 0; i < integrators.Size(); i++)
+   {
+      integrators[i]->AddMultTransposePA(x, y);
+   }
+}
+
+void SumIntegrator::AssembleMF(const FiniteElementSpace &fes)
+{
+   for (int i = 0; i < integrators.Size(); i++)
+   {
+      integrators[i]->AssembleMF(fes);
+   }
+}
+
+void SumIntegrator::AddMultMF(const Vector& x, Vector& y) const
+{
+   for (int i = 0; i < integrators.Size(); i++)
+   {
+      integrators[i]->AddMultTransposeMF(x, y);
+   }
+}
+
+void SumIntegrator::AddMultTransposeMF(const Vector &x, Vector &y) const
+{
+   for (int i = 0; i < integrators.Size(); i++)
+   {
+      integrators[i]->AddMultMF(x, y);
+   }
+}
+
+void SumIntegrator::AssembleDiagonalMF(Vector &diag)
+{
+   for (int i = 0; i < integrators.Size(); i++)
+   {
+      integrators[i]->AssembleDiagonalMF(diag);
+   }
+}
+
+void SumIntegrator::AssembleEA(const FiniteElementSpace &fes, Vector &emat,
+                               const bool add)
+{
+   for (int i = 0; i < integrators.Size(); i++)
+   {
+      integrators[i]->AssembleEA(fes, emat, add);
+   }
+}
+
+void SumIntegrator::AssembleEAInteriorFaces(const FiniteElementSpace &fes,
+                                            Vector &ea_data_int,
+                                            Vector &ea_data_ext,
+                                            const bool add)
+{
+   for (int i = 0; i < integrators.Size(); i++)
+   {
+      integrators[i]->AssembleEAInteriorFaces(fes,ea_data_int,ea_data_ext,add);
+   }
+}
+
+void SumIntegrator::AssembleEABoundaryFaces(const FiniteElementSpace &fes,
+                                            Vector &ea_data_bdr,
+                                            const bool add)
+{
+   for (int i = 0; i < integrators.Size(); i++)
+   {
+      integrators[i]->AssembleEABoundaryFaces(fes, ea_data_bdr, add);
+   }
+}
+
 SumIntegrator::~SumIntegrator()
 {
    if (own_integrators)

fem/bilininteg.hpp

@@ -355,7 +355,7 @@ class SumIntegrator : public BilinearFormIntegrator
 {
 private:
    int own_integrators;
-   DenseMatrix elem_mat;
+   mutable DenseMatrix elem_mat;
    Array<BilinearFormIntegrator*> integrators;
 
 public:
@@ -367,6 +367,55 @@ public:
    virtual void AssembleElementMatrix(const FiniteElement &el,
                                       ElementTransformation &Trans,
                                       DenseMatrix &elmat);
+   virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
+                                       const FiniteElement &test_fe,
+                                       ElementTransformation &Trans,
+                                       DenseMatrix &elmat);
+
+   using BilinearFormIntegrator::AssembleFaceMatrix;
+   virtual void AssembleFaceMatrix(const FiniteElement &el1,
+                                   const FiniteElement &el2,
+                                   FaceElementTransformations &Trans,
+                                   DenseMatrix &elmat);
+
+   virtual void AssembleFaceMatrix(const FiniteElement &trial_face_fe,
+                                   const FiniteElement &test_fe1,
+                                   const FiniteElement &test_fe2,
+                                   FaceElementTransformations &Trans,
+                                   DenseMatrix &elmat);
+
+   using BilinearFormIntegrator::AssemblePA;
+   virtual void AssemblePA(const FiniteElementSpace& fes);
+
+   virtual void AssembleDiagonalPA(Vector &diag);
+
+   virtual void AssemblePAInteriorFaces(const FiniteElementSpace &fes);
+
+   virtual void AssemblePABoundaryFaces(const FiniteElementSpace &fes);
+
+   virtual void AddMultTransposePA(const Vector &x, Vector &y) const;
+
+   virtual void AddMultPA(const Vector& x, Vector& y) const;
+
+   virtual void AssembleMF(const FiniteElementSpace &fes);
+
+   virtual void AddMultMF(const Vector &x, Vector &y) const;
+
+   virtual void AddMultTransposeMF(const Vector &x, Vector &y) const;
+
+   virtual void AssembleDiagonalMF(Vector &diag);
+
+   virtual void AssembleEA(const FiniteElementSpace &fes, Vector &emat,
+                           const bool add);
+
+   virtual void AssembleEAInteriorFaces(const FiniteElementSpace &fes,
+                                        Vector &ea_data_int,
+                                        Vector &ea_data_ext,
+                                        const bool add);
+
+   virtual void AssembleEABoundaryFaces(const FiniteElementSpace &fes,
+                                        Vector &ea_data_bdr,
+                                        const bool add);
 
    virtual ~SumIntegrator();
 };
@@ -1996,6 +2045,8 @@ public:
 
    virtual void AddMultPA(const Vector&, Vector&) const;
 
+   virtual void AddMultTransposePA(const Vector&, Vector&) const;
+
    static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
                                          const FiniteElement &test_fe);
 };
@@ -2059,6 +2110,8 @@ public:
 
    virtual void AddMultPA(const Vector&, Vector&) const;
 
+   virtual void AddMultTransposePA(const Vector&, Vector&) const;
+
    static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
                                          const FiniteElement &test_fe,
                                          ElementTransformation &Trans);
@@ -2120,6 +2173,17 @@ public:
                                          ElementTransformation &Trans);
 };
 
+// Alias for @ConvectionIntegrator.
+using NonconservativeConvectionIntegrator = ConvectionIntegrator;
+
+/// -alpha (u, q . grad v), negative transpose of ConvectionIntegrator
+class ConservativeConvectionIntegrator : public TransposeIntegrator
+{
+public:
+   ConservativeConvectionIntegrator(VectorCoefficient &q, double a = 1.0)
+      : TransposeIntegrator(new ConvectionIntegrator(q, -a)) { }
+};
+
 /// alpha (q . grad u, v) using the "group" FE discretization
 class GroupConvectionIntegrator : public BilinearFormIntegrator
 {
@@ -2684,15 +2748,15 @@ public:
     One use case for this integrator is to discretize the operator -u.grad(v)
     with a DG formulation. The resulting formulation uses the
     ConvectionIntegrator (with coefficient u, and parameter alpha = -1) and the
-    transpose of the DGTraceIntegrator (with coefficient u, and parameters
-    alpha = 1, beta = -1/2 to use the upwind face flux). This discretization and
-    the handling of the inflow and outflow boundaries is illustrated in Example
-    9/9p.
+    transpose of the DGTraceIntegrator (with coefficient u, and parameters alpha
+    = 1, beta = -1/2 to use the upwind face flux, see also
+    NonconservativeDGTraceIntegrator). This discretization and the handling of
+    the inflow and outflow boundaries is illustrated in Example 9/9p.
 
     Another use case for this integrator is to discretize the operator -div(u v)
     with a DG formulation. The resulting formulation is conservative and
-    consists of the transpose of the ConvectionIntegrator (with coefficient u,
-    and parameter alpha = 1) plus the DGTraceIntegrator (with coefficient u, and
+    consists of the ConservativeConvectionIntegrator (with coefficient u, and
+    parameter alpha = -1) plus the DGTraceIntegrator (with coefficient u, and
     parameters alpha = -1, beta = -1/2 to use the upwind face flux).
     */
 class DGTraceIntegrator : public BilinearFormIntegrator
@@ -2711,13 +2775,17 @@ private:
    Vector shape1, shape2;
 
 public:
-   /// Construct integrator with rho = 1.
-   DGTraceIntegrator(VectorCoefficient &_u, double a, double b)
-   { rho = NULL; u = &_u; alpha = a; beta = b; }
+   /// Construct integrator with rho = 1, b = 0.5*a.
+   DGTraceIntegrator(VectorCoefficient &u_, double a)
+   { rho = NULL; u = &u_; alpha = a; beta = 0.5*a; }
 
-   DGTraceIntegrator(Coefficient &_rho, VectorCoefficient &_u,
+   /// Construct integrator with rho = 1.
+   DGTraceIntegrator(VectorCoefficient &u_, double a, double b)
+   { rho = NULL; u = &u_; alpha = a; beta = b; }
+
+   DGTraceIntegrator(Coefficient &_rho, VectorCoefficient &u_,
                      double a, double b)
-   { rho = &_rho; u = &_u; alpha = a; beta = b; }
+   { rho = &_rho; u = &u_; alpha = a; beta = b; }
 
    using BilinearFormIntegrator::AssembleFaceMatrix;
    virtual void AssembleFaceMatrix(const FiniteElement &el1,
@@ -2751,6 +2819,30 @@ private:
    void SetupPA(const FiniteElementSpace &fes, FaceType type);
 };
 
+// Alias for @a DGTraceIntegrator.
+using ConservativeDGTraceIntegrator = DGTraceIntegrator;
+
+/** Integrator that represents the face terms used for the non-conservative
+    DG discretization of the convection equation:
+    -alpha < rho_u (u.n) {v},[w] > + beta < rho_u |u.n| [v],[w] >.
+
+    This integrator can be used with together with ConvectionIntegrator to
+    implement an upwind DG discretization in non-conservative form, see ex9 and
+    ex9p. */
+class NonconservativeDGTraceIntegrator : public TransposeIntegrator
+{
+public:
+   NonconservativeDGTraceIntegrator(VectorCoefficient &u, double a)
+      : TransposeIntegrator(new DGTraceIntegrator(u, -a, 0.5*a)) { }
+
+   NonconservativeDGTraceIntegrator(VectorCoefficient &u, double a, double b)
+      : TransposeIntegrator(new DGTraceIntegrator(u, -a, b)) { }
+
+   NonconservativeDGTraceIntegrator(Coefficient &rho, VectorCoefficient &u,
+                                    double a, double b)
+      : TransposeIntegrator(new DGTraceIntegrator(rho, u, -a, b)) { }
+};
+
 /** Integrator for the DG form:
 
     - < {(Q grad(u)).n}, [v] > + sigma < [u], {(Q grad(v)).n} >

fem/bilininteg_diffusion_pa.cpp

@@ -1904,4 +1904,17 @@ void DiffusionIntegrator::AddMultPA(const Vector &x, Vector &y) const
    }
 }
 
+void DiffusionIntegrator::AddMultTransposePA(const Vector &x, Vector &y) const
+{
+   if (symmetric)
+   {
+      AddMultPA(x, y);
+   }
+   else
+   {
+      MFEM_ABORT("DiffusionIntegrator::AddMultTransposePA only implemented in "
+                 "the symmetric case.")
+   }
+}
+
 } // namespace mfem

fem/bilininteg_mass_pa.cpp

@@ -1226,4 +1226,10 @@ void MassIntegrator::AddMultPA(const Vector &x, Vector &y) const
    }
 }
 
+void MassIntegrator::AddMultTransposePA(const Vector &x, Vector &y) const
+{
+   // Mass integrator is symmetric
+   AddMultPA(x, y);
+}
+
 } // namespace mfem

fem/fe.cpp

@@ -7797,7 +7797,10 @@ NodalTensorFiniteElement::NodalTensorFiniteElement(const int dims,
                                                    const DofMapType dmtype)
    : NodalFiniteElement(dims, GetTensorProductGeometry(dims), Pow(p + 1, dims),
                         p, dims > 1 ? FunctionSpace::Qk : FunctionSpace::Pk),
-     TensorBasisElement(dims, p, VerifyNodal(btype), dmtype) { }
+     TensorBasisElement(dims, p, VerifyNodal(btype), dmtype)
+{
+   lex_ordering = dof_map;
+}
 
 
 PositiveTensorFiniteElement::PositiveTensorFiniteElement(
@@ -8481,23 +8484,33 @@ H1_TriangleElement::H1_TriangleElement(const int p, const int btype)
    Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
 #endif
 
+   int p2p3 = 2*p + 3;
+   auto idx = [p2p3](int i, int j) { return ((p2p3-j)*j)/2+i; };
+   lex_ordering.SetSize(dof);
+
    // vertices
+   lex_ordering[idx(0,0)] = 0;
    Nodes.IntPoint(0).Set2(cp[0], cp[0]);
+   lex_ordering[idx(p,0)] = 1;
    Nodes.IntPoint(1).Set2(cp[p], cp[0]);
+   lex_ordering[idx(0,p)] = 2;
    Nodes.IntPoint(2).Set2(cp[0], cp[p]);
 
    // edges
    int o = 3;
    for (int i = 1; i < p; i++)
    {
+      lex_ordering[idx(i,0)] = o;
       Nodes.IntPoint(o++).Set2(cp[i], cp[0]);
    }
    for (int i = 1; i < p; i++)
    {
+      lex_ordering[idx(p-i,i)] = o;
       Nodes.IntPoint(o++).Set2(cp[p-i], cp[i]);
    }
    for (int i = 1; i < p; i++)
    {
+      lex_ordering[idx(0,p-i)] = o;
       Nodes.IntPoint(o++).Set2(cp[0], cp[p-i]);
    }
 
@@ -8506,6 +8519,7 @@ H1_TriangleElement::H1_TriangleElement(const int p, const int btype)
       for (int i = 1; i + j < p; i++)
       {
          const double w = cp[i] + cp[j] + cp[p-i-j];
+         lex_ordering[idx(i,j)] = o;
          Nodes.IntPoint(o++).Set2(cp[i]/w, cp[j]/w);
       }
 
@@ -8639,36 +8653,56 @@ H1_TetrahedronElement::H1_TetrahedronElement(const int p, const int btype)
    Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
 #endif
 
+   auto tri = [](int k) { return (k*(k + 1))/2; };
+   auto tet = [](int k) { return (k*(k + 1)*(k + 2))/6; };
+   int ndof = tet(p+1);
+   auto idx = [tri, tet, p, ndof](int i, int j, int k)
+   {
+      return ndof - tet(p - k) - tri(p + 1 - k - j) + i;
+   };
+
+   lex_ordering.SetSize(dof);
+
    // vertices
+   lex_ordering[idx(0,0,0)] = 0;
    Nodes.IntPoint(0).Set3(cp[0], cp[0], cp[0]);
+   lex_ordering[idx(p,0,0)] = 1;
    Nodes.IntPoint(1).Set3(cp[p], cp[0], cp[0]);
+   lex_ordering[idx(0,p,0)] = 2;
    Nodes.IntPoint(2).Set3(cp[0], cp[p], cp[0]);
+   lex_ordering[idx(0,0,p)] = 3;
    Nodes.IntPoint(3).Set3(cp[0], cp[0], cp[p]);
 
    // edges (see Tetrahedron::edges in mesh/tetrahedron.cpp)
    int o = 4;
    for (int i = 1; i < p; i++)  // (0,1)
    {
+      lex_ordering[idx(i,0,0)] = o;
       Nodes.IntPoint(o++).Set3(cp[i], cp[0], cp[0]);
    }
    for (int i = 1; i < p; i++)  // (0,2)
    {
+      lex_ordering[idx(0,i,0)] = o;
       Nodes.IntPoint(o++).Set3(cp[0], cp[i], cp[0]);
    }
    for (int i = 1; i < p; i++)  // (0,3)
    {
+      lex_ordering[idx(0,0,i)] = o;
       Nodes.IntPoint(o++).Set3(cp[0], cp[0], cp[i]);
    }
    for (int i = 1; i < p; i++)  // (1,2)
    {
+      lex_ordering[idx(p-i,i,0)] = o;
       Nodes.IntPoint(o++).Set3(cp[p-i], cp[i], cp[0]);
    }
    for (int i = 1; i < p; i++)  // (1,3)
    {
+      lex_ordering[idx(p-i,0,i)] = o;
       Nodes.IntPoint(o++).Set3(cp[p-i], cp[0], cp[i]);
    }
    for (int i = 1; i < p; i++)  // (2,3)
    {
+      lex_ordering[idx(0,p-i,i)] = o;
       Nodes.IntPoint(o++).Set3(cp[0], cp[p-i], cp[i]);
    }
 
@@ -8676,24 +8710,28 @@ H1_TetrahedronElement::H1_TetrahedronElement(const int p, const int btype)
    for (int j = 1; j < p; j++)
       for (int i = 1; i + j < p; i++)  // (1,2,3)
       {
+         lex_ordering[idx(p-i-j,i,j)] = o;
          double w = cp[i] + cp[j] + cp[p-i-j];
          Nodes.IntPoint(o++).Set3(cp[p-i-j]/w, cp[i]/w, cp[j]/w);
       }
    for (int j = 1; j < p; j++)
       for (int i = 1; i + j < p; i++)  // (0,3,2)
       {
+         lex_ordering[idx(0,j,i)] = o;
          double w = cp[i] + cp[j] + cp[p-i-j];
          Nodes.IntPoint(o++).Set3(cp[0], cp[j]/w, cp[i]/w);
       }
    for (int j = 1; j < p; j++)
       for (int i = 1; i + j < p; i++)  // (0,1,3)
       {
+         lex_ordering[idx(i,0,j)] = o;
          double w = cp[i] + cp[j] + cp[p-i-j];
          Nodes.IntPoint(o++).Set3(cp[i]/w, cp[0], cp[j]/w);
       }
    for (int j = 1; j < p; j++)
       for (int i = 1; i + j < p; i++)  // (0,2,1)
       {
+         lex_ordering[idx(j,i,0)] = o;
          double w = cp[i] + cp[j] + cp[p-i-j];
          Nodes.IntPoint(o++).Set3(cp[j]/w, cp[i]/w, cp[0]);
       }
@@ -8703,6 +8741,7 @@ H1_TetrahedronElement::H1_TetrahedronElement(const int p, const int btype)
       for (int j = 1; j + k < p; j++)
          for (int i = 1; i + j + k < p; i++)
          {
+            lex_ordering[idx(i,j,k)] = o;
             double w = cp[i] + cp[j] + cp[k] + cp[p-i-j-k];
             Nodes.IntPoint(o++).Set3(cp[i]/w, cp[j]/w, cp[k]/w);
          }
@@ -9247,7 +9286,22 @@ H1_WedgeElement::H1_WedgeElement(const int p,
    t_dof.SetSize(dof);
    s_dof.SetSize(dof);
 
+   int p2p3 = 2*p + 3, ntri = ((p + 1)*(p + 2))/2;
+   auto idx = [p2p3,ntri](int i, int j, int k)
+   {
+      return k*ntri + ((p2p3-j)*j)/2+i;
+   };
+
+   lex_ordering.SetSize(dof);
+   int o = 0;
+
    // Nodal DoFs
+   lex_ordering[idx(0,0,0)] = o++;
+   lex_ordering[idx(p,0,0)] = o++;
+   lex_ordering[idx(0,p,0)] = o++;
+   lex_ordering[idx(0,0,p)] = o++;
+   lex_ordering[idx(p,0,p)] = o++;
+   lex_ordering[idx(0,p,p)] = o++;
    t_dof[0] = 0; s_dof[0] = 0;
    t_dof[1] = 1; s_dof[1] = 0;
    t_dof[2] = 2; s_dof[2] = 0;
@@ -9256,9 +9310,19 @@ H1_WedgeElement::H1_WedgeElement(const int p,
    t_dof[5] = 2; s_dof[5] = 1;
 
    // Edge DoFs
+   int k = 0;
    int ne = p-1;
    for (int i=1; i<p; i++)
    {
+      lex_ordering[idx(i,0,0)] = o + 0*ne + k;
+      lex_ordering[idx(p-i,i,0)] = o + 1*ne + k;
+      lex_ordering[idx(0,p-i,0)] = o + 2*ne + k;
+      lex_ordering[idx(i,0,p)] = o + 3*ne + k;
+      lex_ordering[idx(p-i,i,p)] = o + 4*ne + k;
+      lex_ordering[idx(0,p-i,p)] = o + 5*ne + k;
+      lex_ordering[idx(0,0,i)] = o + 6*ne + k;
+      lex_ordering[idx(p,0,i)] = o + 7*ne + k;
+      lex_ordering[idx(0,p,i)] = o + 8*ne + k;
       t_dof[5 + 0 * ne + i] = 2 + 0 * ne + i; s_dof[5 + 0 * ne + i] = 0;
       t_dof[5 + 1 * ne + i] = 2 + 1 * ne + i; s_dof[5 + 1 * ne + i] = 0;
       t_dof[5 + 2 * ne + i] = 2 + 2 * ne + i; s_dof[5 + 2 * ne + i] = 0;
@@ -9268,21 +9332,26 @@ H1_WedgeElement::H1_WedgeElement(const int p,
       t_dof[5 + 6 * ne + i] = 0;              s_dof[5 + 6 * ne + i] = i + 1;
       t_dof[5 + 7 * ne + i] = 1;              s_dof[5 + 7 * ne + i] = i + 1;
       t_dof[5 + 8 * ne + i] = 2;              s_dof[5 + 8 * ne + i] = i + 1;
+      ++k;
    }
+   o += 9*ne;
 
    // Triangular Face DoFs
-   int k=0;
+   k=0;
    int nt = (p-1)*(p-2)/2;
    for (int j=1; j<p; j++)
    {
       for (int i=1; i<p-j; i++)
       {
          int l = j - p + (((2 * p - 1) - i) * i) / 2;
+         lex_ordering[idx(i,j,0)] = o+l;
+         lex_ordering[idx(i,j,p)] = o+nt+k;
          t_dof[6 + 9 * ne + k]      = 3 * p + l; s_dof[6 + 9 * ne + k]      = 0;
          t_dof[6 + 9 * ne + nt + k] = 3 * p + k; s_dof[6 + 9 * ne + nt + k] = 1;
          k++;
       }
    }
+   o += 2*nt;
 
    // Quadrilateral Face DoFs
    k=0;
@@ -9291,6 +9360,10 @@ H1_WedgeElement::H1_WedgeElement(const int p,
    {
       for (int i=1; i<p; i++)
       {
+         lex_ordering[idx(i,0,j)] = o+k;
+         lex_ordering[idx(p-i,i,j)] = o+nq+k;
+         lex_ordering[idx(0,p-i,j)] = o+2*nq+k;
+
          t_dof[6 + 9 * ne + 2 * nt + 0 * nq + k] = 2 + 0 * ne + i;
          t_dof[6 + 9 * ne + 2 * nt + 1 * nq + k] = 2 + 1 * ne + i;
          t_dof[6 + 9 * ne + 2 * nt + 2 * nq + k] = 2 + 2 * ne + i;
@@ -9302,6 +9375,7 @@ H1_WedgeElement::H1_WedgeElement(const int p,
          k++;
       }
    }
+   o += 3*nq;
 
    // Interior DoFs
    int m=0;
@@ -9310,8 +9384,9 @@ H1_WedgeElement::H1_WedgeElement(const int p,
       int l=0;
       for (int j=1; j<p; j++)
       {
-         for (int i=1; i<j; i++)
+         for (int i=1; i+j<p; i++)
          {
+            lex_ordering[idx(i,j,k)] = o++;
             t_dof[6 + 9 * ne + 2 * nt + 3 * nq + m] = 3 * p + l;
             s_dof[6 + 9 * ne + 2 * nt + 3 * nq + m] = 1 + k;
             l++; m++;

fem/fe.hpp

@@ -698,6 +698,7 @@ public:
 class NodalFiniteElement : public ScalarFiniteElement
 {
 protected:
+   Array<int> lex_ordering;
    void ProjectCurl_2D(const FiniteElement &fe,
                        ElementTransformation &Trans,
                        DenseMatrix &curl) const;
@@ -746,6 +747,29 @@ public:
    virtual void ProjectDiv(const FiniteElement &fe,
                            ElementTransformation &Trans,
                            DenseMatrix &div) const;
+
+   /** @brief Get an Array<int> that maps lexicographically ordered indices to
+       the indices of the respective nodes/dofs/basis functions. Lexicographic
+       ordering of nodes is defined in terms of reference-space coordinates
+       (x,y,z). Lexicographically ordered nodes are listed first in order of
+       increasing x-coordinate, and then in order of increasing y-coordinate,
+       and finally in order of increasing z-coordinate.
+
+       For example, the six nodes of a quadratic triangle are lexicographically
+       ordered as follows:
+
+       5
+       |\
+       3 4
+       |  \
+       0-1-2
+
+       The resulting array may be empty if the DOFs are already ordered
+       lexicographically, or if the finite element does not support creating
+       this permutation. The array returned is the same as the array given by
+       TensorBasisElement::GetDofMap, but it is also available for non-tensor
+       elements. */
+   const Array<int> &GetLexicographicOrdering() const { return lex_ordering; }
 };
 
 /** @brief Class for finite elements utilizing the

fem/fe_coll.cpp

@@ -1854,19 +1854,16 @@ const int *H1_FECollection::GetDofMap(Geometry::Type GeomType) const
 {
    const int *dof_map = NULL;
    const FiniteElement *fe = H1_Elements[GeomType];
-   switch (GeomType)
+   const NodalFiniteElement *nodal_fe =
+      dynamic_cast<const NodalFiniteElement*>(fe);
+   if (nodal_fe)
    {
-      case Geometry::SEGMENT:
-      case Geometry::SQUARE:
-      case Geometry::CUBE:
-         dof_map = dynamic_cast<const TensorBasisElement *>(fe)
-                   ->GetDofMap().GetData();
-         break;
-      default:
-         MFEM_ABORT("Geometry type " << Geometry::Name[GeomType] << " is not "
-                    "implemented");
-         // The "Cartesian" ordering for other geometries is defined by the
-         // class GeometryRefiner.
+      dof_map = nodal_fe->GetLexicographicOrdering().GetData();
+   }
+   else
+   {
+      MFEM_ABORT("Geometry type " << Geometry::Name[GeomType] << " is not "
+                 "implemented");
    }
    return dof_map;
 }

fem/fem.hpp

@@ -41,6 +41,7 @@
 #include "transfer.hpp"
 #include "fespacehierarchy.hpp"
 #include "multigrid.hpp"
+#include "prmnonlinearform.hpp"
 
 #ifdef MFEM_USE_MPI
 #include "pfespace.hpp"
@@ -48,6 +49,7 @@
 #include "plinearform.hpp"
 #include "pbilinearform.hpp"
 #include "pnonlinearform.hpp"
+#include "pprmnonlinearform.hpp"
 #endif
 
 #ifdef MFEM_USE_SIDRE

fem/geom.cpp

@@ -1165,24 +1165,31 @@ RefinedGeometry * GeometryRefiner::Refine(Geometry::Type Geom,
          Array<int> vi((n+1)*(n+1)*(n+1));
          vi = -1;
          m = 0;
-         for (k = 0; k <= n; k++)
-            for (j = 0; j <= k; j++)
-               for (i = 0; i <= j; i++)
+
+         // vertices are given in lexicographic ordering on the reference
+         // element
+         for (int kk = 0; kk <= n; kk++)
+            for (int jj = 0; jj <= n-kk; jj++)
+               for (int ii = 0; ii <= n-jj-kk; ii++)
                {
                   IntegrationPoint &ip = RG->RefPts.IntPoint(m);
-                  // map the coordinates to the reference tetrahedron
-                  // (0,0,0) -> (0,0,0)
-                  // (0,0,1) -> (1,0,0)
-                  // (1,1,1) -> (0,1,0)
-                  // (0,1,1) -> (0,0,1)
-                  double w = cp[k-j] + cp[i] + cp[j-i] + cp[Times-k];
-                  ip.x = cp[k-j]/w;
-                  ip.y = cp[i]/w;
-                  ip.z = cp[j-i]/w;
+                  double w = cp[ii] + cp[jj] + cp[kk] + cp[Times-ii-jj-kk];
+                  ip.x = cp[ii]/w;
+                  ip.y = cp[jj]/w;
+                  ip.z = cp[kk]/w;
+                  // (ii,jj,kk) are coordinates in the reference tetrahedron,
+                  // transform to coordinates (i,j,k) in the auxiliary
+                  // tetrahedron defined by (0,0,0), (0,0,1), (1,1,1), (0,1,1)
+                  int i = jj;
+                  int j = jj+kk;
+                  int k = ii+jj+kk;
                   l = i + (j + k * (n+1)) * (n+1);
+                  // map from linear Cartesian hex index in the auxiliary tet
+                  // to lexicographic in the reference tet
                   vi[l] = m;
                   m++;
                }
+
          if (m != (n+3)*(n+2)*(n+1)/6)
          {
             mfem_error("GeometryRefiner::Refine() for TETRAHEDRON #1");
@@ -1269,18 +1276,9 @@ RefinedGeometry * GeometryRefiner::Refine(Geometry::Type Geom,
                for (i = 0; i <= n-j; i++, l++)
                {
                   IntegrationPoint &ip = RG->RefPts.IntPoint(l);
-                  if (type == 0)
-                  {
-                     ip.x = double(i) / n;
-                     ip.y = double(j) / n;
-                     ip.z = double(k) / n;
-                  }
-                  else
-                  {
-                     ip.x = cp[i]/(cp[i] + cp[j] + cp[n-i-j]);
-                     ip.y = cp[j]/(cp[i] + cp[j] + cp[n-i-j]);
-                     ip.z = cp[k];
-                  }
+                  ip.x = cp[i]/(cp[i] + cp[j] + cp[n-i-j]);
+                  ip.y = cp[j]/(cp[i] + cp[j] + cp[n-i-j]);
+                  ip.z = cp[k];
                   m++;
                }
          if (m != (n+1)*(n+1)*(n+2)/2)

fem/lininteg.hpp

@@ -406,6 +406,10 @@ private:
    Vector shape;
 
 public:
+   BoundaryFlowIntegrator(Coefficient &_f, VectorCoefficient &_u,
+                          double a)
+   { f = &_f; u = &_u; alpha = a; beta = 0.5*a; }
+
    BoundaryFlowIntegrator(Coefficient &_f, VectorCoefficient &_u,
                           double a, double b)
    { f = &_f; u = &_u; alpha = a; beta = b; }

fem/nonlininteg.cpp

@@ -128,6 +128,92 @@ double BlockNonlinearFormIntegrator::GetElementEnergy(
    return 0.0;
 }
 
+double PrmBlockNonlinearFormIntegrator::GetElementEnergy(const Array<const FiniteElement *> &el,
+                                                         const Array<const FiniteElement *> &pel,
+                                                         ElementTransformation &Tr,
+                                                         const Array<const Vector *> &elfun,
+                                                         const Array<const Vector *> &pelfun)
+{
+    mfem_error("PrmBlockNonlinearFormIntegrator::GetElementEnergy"
+               " is not overloaded!");
+    return 0.0;
+}
+
+void PrmBlockNonlinearFormIntegrator::AssembleFaceGrad(const Array<const FiniteElement *> &el1,
+                                                       const Array<const FiniteElement *> &el2,
+                                                       const Array<const FiniteElement *> &pel1,
+                                                       const Array<const FiniteElement *> &pel2,
+                                                       FaceElementTransformations &Tr,
+                                                       const Array<const Vector *> &elfun,
+                                                       const Array<const Vector *> &pelfun,
+                                                       const Array2D<DenseMatrix *> &elmats)
+{
+    mfem_error("PrmBlockNonlinearFormIntegrator::AssembleFaceGrad"
+               " is not overloaded!");
+}
+
+void PrmBlockNonlinearFormIntegrator::AssembleElementGrad(const Array<const FiniteElement *> &el,
+                                                          const Array<const FiniteElement *> &pel,
+                                                          ElementTransformation &Tr,
+                                                          const Array<const Vector *> &elfun,
+                                                          const Array<const Vector *> &pelfun,
+                                                          const Array2D<DenseMatrix *> &elmats)
+{
+    mfem_error("PrmBlockNonlinearFormIntegrator::AssembleElementGrad"
+               " is not overloaded!");
+}
+
+void PrmBlockNonlinearFormIntegrator::AssembleElementVector(const Array<const FiniteElement *> &el,
+                                                            const Array<const FiniteElement *> &pel,
+                                                            ElementTransformation &Tr,
+                                                            const Array<const Vector *> &elfun,
+                                                            const Array<const Vector *> &pelfun,
+                                                            const Array<Vector *> &elvec)
+{
+    mfem_error("PrmBlockNonlinearFormIntegrator::AssembleElementVector"
+               " is not overloaded!");
+}
+
+void PrmBlockNonlinearFormIntegrator::AssembleFaceVector(const Array<const FiniteElement *> &el1,
+                                                         const Array<const FiniteElement *> &el2,
+                                                         const Array<const FiniteElement *> &pel1,
+                                                         const Array<const FiniteElement *> &pel2,
+                                                         FaceElementTransformations &Tr,
+                                                         const Array<const Vector *> &elfun,
+                                                         const Array<const Vector *> &pelfun,
+                                                         const Array<Vector *> &elvect)
+{
+    mfem_error("PrmBlockNonlinearFormIntegrator::AssembleFaceVector"
+               " is not overloaded!");
+}
+
+
+void PrmBlockNonlinearFormIntegrator::AssemblePrmElementVector(const Array<const FiniteElement *> &el,
+                                                               const Array<const FiniteElement *> &pel,
+                                                               ElementTransformation &Tr,
+                                                               const Array<const Vector *> &elfun,
+                                                               const Array<const Vector *> &alfun,
+                                                               const Array<const Vector *> &pelfun,
+                                                               const Array<Vector *> &elvec)
+{
+    mfem_error("PrmBlockNonlinearFormIntegrator::AssemblePrmElementVector"
+               " is not overloaded!");
+}
+
+void PrmBlockNonlinearFormIntegrator::AssemblePrmFaceVector(const Array<const FiniteElement *> &el1,
+                                                         const Array<const FiniteElement *> &el2,
+                                                         const Array<const FiniteElement *> &pel1,
+                                                         const Array<const FiniteElement *> &pel2,
+                                                         FaceElementTransformations &Tr,
+                                                         const Array<const Vector *> &elfun,
+                                                         const Array<const Vector *> &alfun,
+                                                         const Array<const Vector *> &pelfun,
+                                                         const Array<Vector *> &elvect)
+{
+    mfem_error("PrmBlockNonlinearFormIntegrator::AssemblePrmFaceVector"
+               " is not overloaded!");
+}
+
 
 double InverseHarmonicModel::EvalW(const DenseMatrix &J) const
 {

fem/nonlininteg.hpp

@@ -130,6 +130,80 @@ public:
 };
 
 
+/** The abstract base class PrmBlockNonlinearFormIntegrator is
+    a generalization of the BlockNonlinearFormIntegrator class suitable
+    for block state and parameter vectors. */
+class PrmBlockNonlinearFormIntegrator
+{
+public:
+   /// Compute the local energy
+   virtual double GetElementEnergy(const Array<const FiniteElement *>&el,
+                                   const Array<const FiniteElement *>&pel,
+                                   ElementTransformation &Tr,
+                                   const Array<const Vector *>&elfun,
+                                   const Array<const Vector *>&pelfun);
+
+   /// Perform the local action of the BlockNonlinearFormIntegrator
+   virtual void AssembleElementVector(const Array<const FiniteElement *> &el,
+                                      const Array<const FiniteElement *>&pel,
+                                      ElementTransformation &Tr,
+                                      const Array<const Vector *> &elfun,
+                                      const Array<const Vector *>&pelfun,
+                                      const Array<Vector *> &elvec);
+
+   virtual void AssembleFaceVector(const Array<const FiniteElement *> &el1,
+                                   const Array<const FiniteElement *> &el2,
+                                   const Array<const FiniteElement *> &pel1,
+                                   const Array<const FiniteElement *> &pel2,
+                                   FaceElementTransformations &Tr,
+                                   const Array<const Vector *> &elfun,
+                                   const Array<const Vector *>&pelfun,
+                                   const Array<Vector *> &elvect);
+
+    /// Perform the local action on the parameters of the BlockNonlinearFormIntegrator
+    virtual void AssemblePrmElementVector(const Array<const FiniteElement *> &el,
+                                       const Array<const FiniteElement *>&pel,
+                                       ElementTransformation &Tr,
+                                       const Array<const Vector *> &elfun,
+                                       const Array<const Vector *> &alfun,
+                                       const Array<const Vector *>&pelfun,
+                                       const Array<Vector *> &elvec);
+
+    virtual void AssemblePrmFaceVector(const Array<const FiniteElement *> &el1,
+                                    const Array<const FiniteElement *> &el2,
+                                    const Array<const FiniteElement *> &pel1,
+                                    const Array<const FiniteElement *> &pel2,
+                                    FaceElementTransformations &Tr,
+                                    const Array<const Vector *> &elfun,
+                                    const Array<const Vector *> &alfun,
+                                    const Array<const Vector *>&pelfun,
+                                    const Array<Vector *> &elvect);
+
+
+
+   /// Assemble the local gradient matrix
+   virtual void AssembleElementGrad(const Array<const FiniteElement*> &el,
+                                    const Array<const FiniteElement *>&pel,
+                                    ElementTransformation &Tr,
+                                    const Array<const Vector *> &elfun,
+                                    const Array<const Vector *>&pelfun,
+                                    const Array2D<DenseMatrix *> &elmats);
+
+   virtual void AssembleFaceGrad(const Array<const FiniteElement *>&el1,
+                                 const Array<const FiniteElement *>&el2,
+                                 const Array<const FiniteElement *> &pel1,
+                                 const Array<const FiniteElement *> &pel2,
+                                 FaceElementTransformations &Tr,
+                                 const Array<const Vector *> &elfun,
+                                 const Array<const Vector *>&pelfun,
+                                 const Array2D<DenseMatrix *> &elmats);
+
+
+   virtual ~PrmBlockNonlinearFormIntegrator() { }
+};
+
+
+
 /// Abstract class for hyperelastic models
 class HyperelasticModel
 {

fem/pbilinearform.hpp

@@ -249,6 +249,9 @@ public:
        @a A. */
    void ParallelAssemble(OperatorHandle &A);
 
+   using MixedBilinearForm::FormRectangularSystemMatrix;
+   using MixedBilinearForm::FormRectangularLinearSystem;
+
    /** @brief Return in @a A a parallel (on truedofs) version of this operator.
 
        This returns the same operator as FormRectangularLinearSystem(), but does
@@ -313,6 +316,8 @@ public:
        domain and range parallel finite element spaces */
    void GetParBlocks(Array2D<HypreParMatrix *> &blocks) const;
 
+   using MixedBilinearForm::FormRectangularSystemMatrix;
+
    /** @brief Return in @a A a parallel (on truedofs) version of this operator. */
    virtual void FormRectangularSystemMatrix(OperatorHandle &A);
 

fem/pprmnonlinearform.cpp

@@ -0,0 +1,363 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_MPI
+
+#include "fem.hpp"
+
+
+namespace mfem
+{
+
+ParPrmBlockNonlinearForm::ParPrmBlockNonlinearForm(Array<ParFiniteElementSpace *> &pf,
+                                                   Array<ParFiniteElementSpace *> &ppf)
+    :PrmBlockNonlinearForm()
+{
+    pBlockGrad = nullptr;
+    SetParSpaces(pf,ppf);
+}
+
+void ParPrmBlockNonlinearForm::SetParSpaces(Array<ParFiniteElementSpace *> &pf,
+                                            Array<ParFiniteElementSpace *> &pprmf)
+{
+   delete pBlockGrad;
+   pBlockGrad = nullptr;
+
+   for (int s1=0; s1<fes.Size(); ++s1)
+   {
+      for (int s2=0; s2<fes.Size(); ++s2)
+      {
+         delete phBlockGrad(s1,s2);
+      }
+   }
+
+   Array<FiniteElementSpace *> serialSpaces(pf.Size());
+   Array<FiniteElementSpace *> prmserialSpaces(pprmf.Size());
+   for (int s=0; s<pf.Size(); s++)
+   {
+      serialSpaces[s] = (FiniteElementSpace *) pf[s];
+   }
+   for (int s=0; s<pprmf.Size(); s++)
+   {
+      prmserialSpaces[s] = (FiniteElementSpace *) pprmf[s];
+   }
+
+   SetSpaces(serialSpaces,prmserialSpaces);
+
+   phBlockGrad.SetSize(fes.Size(), fes.Size());
+
+   for (int s1=0; s1<fes.Size(); ++s1)
+   {
+      for (int s2=0; s2<fes.Size(); ++s2)
+      {
+         phBlockGrad(s1,s2) = new OperatorHandle(Operator::Hypre_ParCSR);
+      }
+   }
+}
+
+ParFiniteElementSpace * ParPrmBlockNonlinearForm::ParFESpace(int k)
+{
+   return (ParFiniteElementSpace *)fes[k];
+}
+
+const ParFiniteElementSpace *ParPrmBlockNonlinearForm::ParFESpace(int k) const
+{
+   return (const ParFiniteElementSpace *)fes[k];
+}
+
+
+ParFiniteElementSpace * ParPrmBlockNonlinearForm::ParPrmFESpace(int k)
+{
+   return (ParFiniteElementSpace *)prmfes[k];
+}
+
+const ParFiniteElementSpace *ParPrmBlockNonlinearForm::ParPrmFESpace(int k) const
+{
+   return (const ParFiniteElementSpace *)prmfes[k];
+}
+
+// Here, rhs is a true dof vector
+void ParPrmBlockNonlinearForm::SetEssentialBC(const
+                                           Array<Array<int> *>&bdr_attr_is_ess,
+                                           Array<Vector *> &rhs)
+{
+   Array<Vector *> nullarray(fes.Size());
+   nullarray = NULL;
+
+   PrmBlockNonlinearForm::SetEssentialBC(bdr_attr_is_ess, nullarray);
+
+   for (int s = 0; s < fes.Size(); ++s)
+   {
+      if (rhs[s])
+      {
+         rhs[s]->SetSubVector(*ess_tdofs[s], 0.0);
+      }
+   }
+}
+
+void ParPrmBlockNonlinearForm::SetPrmEssentialBC(const
+                                           Array<Array<int> *>&bdr_attr_is_ess,
+                                           Array<Vector *> &rhs)
+{
+    Array<Vector *> nullarray(fes.Size());
+    nullarray = NULL;
+
+    PrmBlockNonlinearForm::SetPrmEssentialBC(bdr_attr_is_ess, nullarray);
+
+    for (int s = 0; s < prmfes.Size(); ++s)
+    {
+       if (rhs[s])
+       {
+          rhs[s]->SetSubVector(*prmess_tdofs[s], 0.0);
+       }
+    }
+}
+
+double ParPrmBlockNonlinearForm::GetEnergy(const Vector &x) const
+{
+   xs_true.Update(x.GetData(), block_trueOffsets);
+   xs.Update(block_offsets);
+
+   for (int s = 0; s < fes.Size(); ++s)
+   {
+      fes[s]->GetProlongationMatrix()->Mult(xs_true.GetBlock(s), xs.GetBlock(s));
+   }
+
+   double enloc = PrmBlockNonlinearForm::GetEnergyBlocked(xs,xdv);
+   double englo = 0.0;
+
+   MPI_Allreduce(&enloc, &englo, 1, MPI_DOUBLE, MPI_SUM,
+                 ParFESpace(0)->GetComm());
+
+   return englo;
+}
+
+void ParPrmBlockNonlinearForm::Mult(const Vector &x, Vector &y) const
+{
+   xs_true.Update(x.GetData(), block_trueOffsets);
+   ys_true.Update(y.GetData(), block_trueOffsets);
+   xs.Update(block_offsets);
+   ys.Update(block_offsets);
+
+   for (int s=0; s<fes.Size(); ++s)
+   {
+      fes[s]->GetProlongationMatrix()->Mult(
+         xs_true.GetBlock(s), xs.GetBlock(s));
+   }
+
+   PrmBlockNonlinearForm::MultBlocked(xs, xdv, ys);
+
+   if (fnfi.Size() > 0)
+   {
+      MFEM_ABORT("TODO: assemble contributions from shared face terms");
+   }
+
+   for (int s=0; s<fes.Size(); ++s)
+   {
+      fes[s]->GetProlongationMatrix()->MultTranspose(
+         ys.GetBlock(s), ys_true.GetBlock(s));
+
+      ys_true.GetBlock(s).SetSubVector(*ess_tdofs[s], 0.0);
+   }
+}
+
+/// Block T-Vector to Block T-Vector
+void ParPrmBlockNonlinearForm::PrmMult(const Vector &x, Vector &y) const
+{
+    xs_true.Update(x.GetData(), prmblock_trueOffsets);
+    ys_true.Update(y.GetData(), prmblock_trueOffsets);
+    prmxs.Update(prmblock_offsets);
+    prmys.Update(prmblock_offsets);
+
+    for (int s=0; s<prmfes.Size(); ++s)
+    {
+       prmfes[s]->GetProlongationMatrix()->Mult(
+          xs_true.GetBlock(s), prmxs.GetBlock(s));
+    }
+
+    PrmBlockNonlinearForm::MultPrmBlocked(xsv,adv,xdv,prmys);
+
+    if (fnfi.Size() > 0)
+    {
+       MFEM_ABORT("TODO: assemble contributions from shared face terms");
+    }
+
+    for (int s=0; s<prmfes.Size(); ++s)
+    {
+       prmfes[s]->GetProlongationMatrix()->MultTranspose(
+          prmys.GetBlock(s), ys_true.GetBlock(s));
+
+       ys_true.GetBlock(s).SetSubVector(*prmess_tdofs[s], 0.0);
+    }
+
+}
+
+/// Return the local gradient matrix for the given true-dof vector x
+const BlockOperator & ParPrmBlockNonlinearForm::GetLocalGradient(
+   const Vector &x) const
+{
+   xs_true.Update(x.GetData(), block_trueOffsets);
+   xs.Update(block_offsets);
+
+   for (int s=0; s<fes.Size(); ++s)
+   {
+      fes[s]->GetProlongationMatrix()->Mult(
+         xs_true.GetBlock(s), xs.GetBlock(s));
+   }
+
+   PrmBlockNonlinearForm::ComputeGradientBlocked(xs,xdv); // (re)assemble Grad with b.c.
+
+   delete BlockGrad;
+   BlockGrad = new BlockOperator(block_offsets);
+
+   for (int i = 0; i < fes.Size(); ++i)
+   {
+      for (int j = 0; j < fes.Size(); ++j)
+      {
+         BlockGrad->SetBlock(i, j, Grads(i, j));
+      }
+   }
+   return *BlockGrad;
+}
+
+// Set the operator type id for the parallel gradient matrix/operator.
+void ParPrmBlockNonlinearForm::SetGradientType(Operator::Type tid)
+{
+   for (int s1=0; s1<fes.Size(); ++s1)
+   {
+      for (int s2=0; s2<fes.Size(); ++s2)
+      {
+         phBlockGrad(s1,s2)->SetType(tid);
+      }
+   }
+}
+
+BlockOperator & ParPrmBlockNonlinearForm::GetGradient(const Vector &x) const
+{
+   if (pBlockGrad == NULL)
+   {
+      pBlockGrad = new BlockOperator(block_trueOffsets);
+   }
+
+   Array<const ParFiniteElementSpace *> pfes(fes.Size());
+
+   for (int s1=0; s1<fes.Size(); ++s1)
+   {
+      pfes[s1] = ParFESpace(s1);
+
+      for (int s2=0; s2<fes.Size(); ++s2)
+      {
+         phBlockGrad(s1,s2)->Clear();
+      }
+   }
+
+   GetLocalGradient(x); // gradients are stored in 'Grads'
+
+   if (fnfi.Size() > 0)
+   {
+      MFEM_ABORT("TODO: assemble contributions from shared face terms");
+   }
+
+   for (int s1=0; s1<fes.Size(); ++s1)
+   {
+      for (int s2=0; s2<fes.Size(); ++s2)
+      {
+         OperatorHandle dA(phBlockGrad(s1,s2)->Type()),
+                        Ph(phBlockGrad(s1,s2)->Type()),
+                        Rh(phBlockGrad(s1,s2)->Type());
+
+         if (s1 == s2)
+         {
+            dA.MakeSquareBlockDiag(pfes[s1]->GetComm(), pfes[s1]->GlobalVSize(),
+                                   pfes[s1]->GetDofOffsets(), Grads(s1,s1));
+            Ph.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
+            phBlockGrad(s1,s1)->MakePtAP(dA, Ph);
+
+            OperatorHandle Ae;
+            Ae.EliminateRowsCols(*phBlockGrad(s1,s1), *ess_tdofs[s1]);
+         }
+         else
+         {
+            dA.MakeRectangularBlockDiag(pfes[s1]->GetComm(),
+                                        pfes[s1]->GlobalVSize(),
+                                        pfes[s2]->GlobalVSize(),
+                                        pfes[s1]->GetDofOffsets(),
+                                        pfes[s2]->GetDofOffsets(),
+                                        Grads(s1,s2));
+            Rh.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
+            Ph.ConvertFrom(pfes[s2]->Dof_TrueDof_Matrix());
+
+            phBlockGrad(s1,s2)->MakeRAP(Rh, dA, Ph);
+
+            phBlockGrad(s1,s2)->EliminateRows(*ess_tdofs[s1]);
+            phBlockGrad(s1,s2)->EliminateCols(*ess_tdofs[s2]);
+         }
+
+         pBlockGrad->SetBlock(s1, s2, phBlockGrad(s1,s2)->Ptr());
+      }
+   }
+
+   return *pBlockGrad;
+}
+
+ParPrmBlockNonlinearForm::~ParPrmBlockNonlinearForm()
+{
+   delete pBlockGrad;
+   for (int s1=0; s1<fes.Size(); ++s1)
+   {
+      for (int s2=0; s2<fes.Size(); ++s2)
+      {
+         delete phBlockGrad(s1,s2);
+      }
+   }
+}
+
+
+void ParPrmBlockNonlinearForm::SetStateFields(const Vector &xv) const
+{
+   xs_true.Update(xv.GetData(), block_trueOffsets);
+   xsv.Update(block_offsets);
+   for (int s=0; s<fes.Size(); ++s)
+   {
+      fes[s]->GetProlongationMatrix()->Mult(
+         xs_true.GetBlock(s), xsv.GetBlock(s));
+   }
+}
+
+
+void ParPrmBlockNonlinearForm::SetAdjointFields(const Vector &av) const
+{
+    xs_true.Update(av.GetData(), block_trueOffsets);
+    adv.Update(block_offsets);
+    for (int s=0; s<fes.Size(); ++s)
+    {
+       fes[s]->GetProlongationMatrix()->Mult(
+          xs_true.GetBlock(s), adv.GetBlock(s));
+    }
+}
+
+void ParPrmBlockNonlinearForm::SetPrmFields(const Vector &dv) const
+{
+    xs_true.Update(dv.GetData(),prmblock_trueOffsets);
+    xdv.Update(prmblock_offsets);
+    for (int s=0; s<prmfes.Size(); ++s)
+    {
+       prmfes[s]->GetProlongationMatrix()->Mult(
+          xs_true.GetBlock(s), xdv.GetBlock(s));
+    }
+}
+
+
+}
+
+#endif

fem/pprmnonlinearform.hpp

@@ -0,0 +1,105 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_PPRMNONLINEARFORM
+#define MFEM_PPRMNONLINEARFORM
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_MPI
+
+#include "pgridfunc.hpp"
+#include "prmnonlinearform.hpp"
+
+namespace mfem
+{
+
+/** @brief A class representing a general parametric parallel block nonlinear operator
+    defined on the Cartesian product of multiple ParFiniteElementSpace%s. */
+/** The ParPrmBlockNonlinearForm takes as input, and returns as output, vectors on
+    the true dofs. */
+class ParPrmBlockNonlinearForm : public PrmBlockNonlinearForm
+{
+protected:
+   mutable BlockVector xs_true, ys_true;
+   mutable Array2D<OperatorHandle *> phBlockGrad;
+   mutable BlockOperator *pBlockGrad;
+
+public:
+   /// Computes the energy of the system
+   virtual double GetEnergy(const Vector &x) const;
+
+   /// Construct an empty ParPrmBlockNonlinearForm. Initialize with SetParSpaces().
+   ParPrmBlockNonlinearForm() : pBlockGrad(nullptr) { }
+
+   /** @brief Construct a ParPrmBlockNonlinearForm on the given set of
+       parametric and state ParFiniteElementSpace%s. */
+   ParPrmBlockNonlinearForm(Array<ParFiniteElementSpace *> &pf, Array<ParFiniteElementSpace *> &ppf );
+
+   /// Return the @a k-th parallel FE state space of the ParPrmBlockNonlinearForm.
+   ParFiniteElementSpace *ParFESpace(int k);
+   /** @brief Return the @a k-th parallel FE state space of the ParPrmBlockNonlinearForm
+       (const version). */
+   const ParFiniteElementSpace *ParFESpace(int k) const;
+
+   /// Return the @a k-th parallel FE parameters space of the ParPrmBlockNonlinearForm.
+   ParFiniteElementSpace *ParPrmFESpace(int k);
+   /** @brief Return the @a k-th parallel FE parameters space of the ParPrmBlockNonlinearForm
+       (const version). */
+   const ParFiniteElementSpace *ParPrmFESpace(int k) const;
+
+
+   /** @brief After a call to SetParSpaces(), the essential b.c. and the
+       gradient-type (if different from the default) must be set again. */
+   void SetParSpaces(Array<ParFiniteElementSpace *> &pf, Array<ParFiniteElementSpace *> &pprmf);
+
+   // Here, rhs is a true dof vector
+   virtual void SetEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
+                                  Array<Vector *> &rhs);
+
+   // Here, rhs is a true dof vector
+   virtual void SetPrmEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
+                                  Array<Vector *> &rhs);
+
+
+   /// Block T-Vector to Block T-Vector
+   virtual void Mult(const Vector &x, Vector &y) const;
+
+   /// Block T-Vector to Block T-Vector
+   virtual void PrmMult(const Vector &x, Vector &y) const;
+
+   /// Return the local block gradient matrix for the given true-dof vector x
+   const BlockOperator &GetLocalGradient(const Vector &x) const;
+
+   virtual BlockOperator &GetGradient(const Vector &x) const;
+
+   /** @brief Set the operator type id for the blocks of the parallel gradient
+       matrix/operator. The default type is Operator::Hypre_ParCSR. */
+   void SetGradientType(Operator::Type tid);
+
+   /// Destructor.
+   virtual ~ParPrmBlockNonlinearForm();
+
+   /// Set the state fields
+   virtual void SetStateFields(const Vector &xv) const;
+
+   /// Set the adjoint fields
+   virtual void SetAdjointFields(const Vector &av) const;
+
+   /// Set the parameters/design fields
+   virtual void SetPrmFields(const Vector &dv) const;
+
+};
+
+}
+
+#endif
+#endif

fem/prmnonlinearform.hpp

@@ -0,0 +1,231 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_PRMNONLINEARFORM
+#define MFEM_PRMNONLINEARFORM
+
+#include "../config/config.hpp"
+#include "nonlininteg.hpp"
+#include "nonlinearform_ext.hpp"
+#include "bilinearform.hpp"
+#include "gridfunc.hpp"
+
+namespace mfem
+{
+
+/** @brief A class representing a general parametric block nonlinear operator
+    defined on the Cartesian product of multiple FiniteElementSpace%s. */
+class PrmBlockNonlinearForm : public Operator
+{
+protected:
+    /// FE spaces on which the form lives.
+    Array<FiniteElementSpace*> fes;
+
+    /// FE spaces for the parametric fields
+    Array<FiniteElementSpace*> prmfes;
+
+    int prmheight;
+    int prmwidth;
+
+
+    /// Set of Domain Integrators to be assembled (added).
+    Array<PrmBlockNonlinearFormIntegrator*> dnfi;
+
+    /// Set of interior face Integrators to be assembled (added).
+    Array<PrmBlockNonlinearFormIntegrator*> fnfi;
+
+    /// Set of Boundary Face Integrators to be assembled (added).
+    Array<PrmBlockNonlinearFormIntegrator*> bfnfi;
+    Array<Array<int>*>           bfnfi_marker;
+
+    /** Auxiliary block-vectors for wrapping input and output vectors or holding
+           GridFunction-like block-vector data (e.g. in parallel). */
+    mutable BlockVector xs, ys;
+    mutable BlockVector prmxs, prmys;
+
+
+    /** Auxiliary block-vectors for holding
+           GridFunction-like block-vector data (e.g. in parallel). */
+    mutable BlockVector xsv;
+
+    /** Auxiliary block-vectors for holding
+           GridFunction-like block-vector data for the parameter fields
+                                                (e.g. in parallel). */
+    mutable BlockVector xdv;
+    /** Auxiliary block-vectors for holding
+           GridFunction-like block-vector data for the adjoint fields
+                                                (e.g. in parallel). */
+    mutable BlockVector adv;
+
+
+    mutable Array2D<SparseMatrix*> Grads, cGrads;
+    mutable BlockOperator *BlockGrad;
+
+    // A list of the offsets
+    Array<int> block_offsets;
+    Array<int> block_trueOffsets;
+    // A list with the offsets for the parametric fields
+    Array<int> prmblock_offsets;
+    Array<int> prmblock_trueOffsets;
+
+    // Array of Arrays of tdofs for each space in 'fes'
+    Array<Array<int> *> ess_tdofs;
+
+    // Array of Arrays of tdofs for each space in 'prmfes'
+    Array<Array<int> *> prmess_tdofs;
+
+    /// Array of pointers to the prolongation matrix of fes, may be NULL
+    Array<const Operator *> P;
+
+    /// Array of pointers to the prolongation matrix of prmfes, may be NULL
+    Array<const Operator *> Pprm;
+
+    /// Array of results of dynamic-casting P to SparseMatrix pointer
+    Array<const SparseMatrix *> cP;
+
+    /// Array of results of dynamic-casting Pprm to SparseMatrix pointer
+    Array<const SparseMatrix *> cPprm;
+
+
+    /// Indicator if the Operator is part of a parallel run
+    bool is_serial = true;
+
+    /// Indicator if the Operator needs prolongation on assembly
+    bool needs_prolongation = false;
+
+    /// Indicator if the Operator needs prolongation on assembly
+    bool prmneeds_prolongation = false;
+
+
+    mutable BlockVector aux1, aux2;
+
+    mutable BlockVector prmaux1, prmaux2;
+
+    const BlockVector &Prolongate(const BlockVector &bx) const;
+
+    const BlockVector &PrmProlongate(const BlockVector &bx) const;
+
+    /// Specialized version of GetEnergy() for BlockVectors
+    //double GetEnergyBlocked(const BlockVector &bx) const;
+    double GetEnergyBlocked(const BlockVector &bx, const BlockVector &dx) const;
+
+
+    /// Specialized version of Mult() for BlockVector%s
+    /// Block L-Vector to Block L-Vector
+    void MultBlocked(const BlockVector &bx, const BlockVector &dx, BlockVector &by) const;
+
+    /// Specialized version of Mult() for BlockVector%s
+    /// Block L-Vector to Block L-Vector
+    /// bx - state vector, ax - adjoint vector, dx - parametric fields
+    /// dy = ax' d(residual(bx))/d(dx)
+    void MultPrmBlocked(const BlockVector &bx, const BlockVector & ax, const BlockVector &dx, BlockVector &dy) const;
+
+
+    /// Specialized version of GetGradient() for BlockVector
+    //void ComputeGradientBlocked(const BlockVector &bx) const;
+    void ComputeGradientBlocked(const BlockVector &bx, const BlockVector &dx) const;
+
+public:
+   /// Construct an empty BlockNonlinearForm. Initialize with SetSpaces().
+   PrmBlockNonlinearForm();
+
+   /// Construct a BlockNonlinearForm on the given set of FiniteElementSpace%s.
+   PrmBlockNonlinearForm(Array<FiniteElementSpace *> &f, Array<FiniteElementSpace *> &pf );
+
+   /// Return the @a k-th FE space of the PrmBlockNonlinearForm.
+   FiniteElementSpace *FESpace(int k) { return fes[k]; }
+
+   /// Return the @a k-th parametric FE space of the PrmBlockNonlinearForm.
+   FiniteElementSpace *PrmFESpace(int k) { return prmfes[k]; }
+
+
+   /// Return the @a k-th FE space of the BlockNonlinearForm (const version).
+   const FiniteElementSpace *FESpace(int k) const { return fes[k]; }
+
+   /// Return the @a k-th parametric FE space of the BlockNonlinearForm (const version).
+   const FiniteElementSpace *PrmFESpace(int k) const { return prmfes[k]; }
+
+   Array<PrmBlockNonlinearFormIntegrator*>& GetDNFI(){ return dnfi;}
+
+
+   /// (Re)initialize the PrmBlockNonlinearForm.
+   /** After a call to SetSpaces(), the essential b.c. must be set again. */
+   void SetSpaces(Array<FiniteElementSpace *> &f, Array<FiniteElementSpace *> &prmf);
+
+   /// Return the regular dof offsets.
+   const Array<int> &GetBlockOffsets() const { return block_offsets; }
+   /// Return the true-dof offsets.
+   const Array<int> &GetBlockTrueOffsets() const { return block_trueOffsets; }
+
+   /// Return the regular dof offsets for the parameters.
+   const Array<int> &PrmGetBlockOffsets() const { return prmblock_offsets; }
+   /// Return the true-dof offsets for the parameters.
+   const Array<int> &PrmGetBlockTrueOffsets() const { return prmblock_trueOffsets; }
+
+   /// Adds new Domain Integrator.
+   void AddDomainIntegrator(PrmBlockNonlinearFormIntegrator *nlfi)
+   { dnfi.Append(nlfi); }
+
+   /// Adds new Interior Face Integrator.
+   void AddInteriorFaceIntegrator(PrmBlockNonlinearFormIntegrator *nlfi)
+   { fnfi.Append(nlfi); }
+
+   /// Adds new Boundary Face Integrator.
+   void AddBdrFaceIntegrator(PrmBlockNonlinearFormIntegrator *nlfi)
+   { bfnfi.Append(nlfi); bfnfi_marker.Append(NULL); }
+
+   /** @brief Adds new Boundary Face Integrator, restricted to specific boundary
+       attributes. */
+   void AddBdrFaceIntegrator(PrmBlockNonlinearFormIntegrator *nlfi,
+                             Array<int> &bdr_marker);
+
+   virtual void SetEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
+                               Array<Vector *> &rhs);
+
+   virtual void SetPrmEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
+                               Array<Vector *> &rhs);
+
+
+   virtual double GetEnergy(const Vector &x) const;
+
+   /// Method is only called in serial, the parallel version calls MultBlocked
+   /// directly.
+   virtual void Mult(const Vector &x, Vector &y) const;
+
+   /// Method is only called in serial, the parallel version calls MultBlocked
+   /// directly.
+   virtual void PrmMult(const Vector &x, Vector &t) const;
+
+   /// Method is only called in serial, the parallel version calls
+   /// GetGradientBlocked directly.
+   virtual Operator &GetGradient(const Vector &x) const;
+
+   /// Set the state fields
+   virtual void SetStateFields(const Vector &xv) const;
+
+   /// Set the adjoint fields
+   virtual void SetAdjointFields(const Vector &av) const;
+
+   /// Set the parameters/design fields
+   virtual void SetPrmFields(const Vector &dv) const;
+
+
+
+   /// Destructor.
+   virtual ~PrmBlockNonlinearForm();
+
+};
+
+}
+
+#endif
+
+

general/CMakeLists.txt

@@ -27,6 +27,7 @@ list(APPEND SRCS
   stable3d.cpp
   table.cpp
   tic_toc.cpp
+  tinyxml2.cpp
   version.cpp
   )
 
@@ -55,6 +56,7 @@ list(APPEND HDRS
   table.hpp
   tassign.hpp
   tic_toc.hpp
+  tinyxml2.h
   text.hpp
   version.hpp
   )

linalg/hypre.cpp

@@ -128,6 +128,15 @@ HypreParVector::HypreParVector(ParFiniteElementSpace *pfes)
    own_ParVector = 1;
 }
 
+void HypreParVector::WrapHypreParVector(hypre_ParVector *y, bool owner)
+{
+   if (own_ParVector) { hypre_ParVectorDestroy(x); }
+   Destroy();
+   x = y;
+   _SetDataAndSize_();
+   own_ParVector = owner;
+}
+
 Vector * HypreParVector::GlobalVector() const
 {
    hypre_Vector *hv = hypre_ParVectorToVectorAll(*this);
@@ -926,6 +935,15 @@ void HypreParMatrix::GetOffd(SparseMatrix &offd, HYPRE_Int* &cmap) const
    cmap = A->col_map_offd;
 }
 
+void HypreParMatrix::MergeDiagAndOffd(SparseMatrix &merged)
+{
+   SparseMatrix tmp_wrapper;
+   hypre_CSRMatrix *hypre_merged = hypre_MergeDiagAndOffd(A);
+   MakeWrapper(hypre_merged, tmp_wrapper);
+   merged = tmp_wrapper;
+   hypre_CSRMatrixDestroy(hypre_merged);
+}
+
 void HypreParMatrix::GetBlocks(Array2D<HypreParMatrix*> &blocks,
                                bool interleaved_rows,
                                bool interleaved_cols) const
@@ -966,6 +984,46 @@ HypreParMatrix * HypreParMatrix::Transpose() const
    return new HypreParMatrix(At);
 }
 
+#if MFEM_HYPRE_VERSION >= 21800
+HypreParMatrix *HypreParMatrix::ExtractSubmatrix(const Array<int> &indices,
+                                                 double threshhold) const
+{
+   if (!(A->comm))
+   {
+      hypre_MatvecCommPkgCreate(A);
+   }
+
+   hypre_ParCSRMatrix *submat;
+
+   // Get number of rows stored on this processor
+   int local_num_vars = hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(A));
+
+   // Form hypre CF-splitting array designating submatrix as F-points (-1)
+   Array<int> CF_marker(local_num_vars);
+   CF_marker = 1;
+   for (int j=0; j<indices.Size(); j++)
+   {
+      if (indices[j] > local_num_vars)
+      {
+         MFEM_WARNING("WARNING : " << indices[j] << " > " << local_num_vars);
+      }
+      CF_marker[indices[j]] = -1;
+   }
+
+   // Construct cpts_global array on hypre matrix structure
+   HYPRE_BigInt *cpts_global;
+   hypre_BoomerAMGCoarseParms(MPI_COMM_WORLD, local_num_vars, 1, NULL,
+                              CF_marker, NULL, &cpts_global);
+
+   // Extract submatrix into *submat
+   hypre_ParCSRMatrixExtractSubmatrixFC(A, CF_marker, cpts_global,
+                                        "FF", &submat, threshhold);
+
+   mfem_hypre_TFree(cpts_global);
+   return new HypreParMatrix(submat);
+}
+#endif
+
 HYPRE_Int HypreParMatrix::Mult(HypreParVector &x, HypreParVector &y,
                                double a, double b)
 {
@@ -1569,6 +1627,36 @@ void HypreParMatrix::Destroy()
    }
 }
 
+#if MFEM_HYPRE_VERSION >= 21800
+
+void BlockInverseScale(const HypreParMatrix *A, HypreParMatrix *C,
+                       const Vector *b, HypreParVector *d,
+                       int blocksize, BlockInverseScaleJob job)
+{
+   if (job == BlockInverseScaleJob::MATRIX_ONLY ||
+       job == BlockInverseScaleJob::MATRIX_AND_RHS)
+   {
+      hypre_ParCSRMatrix *C_hypre;
+      hypre_ParcsrBdiagInvScal(*A, blocksize, &C_hypre);
+      hypre_ParCSRMatrixDropSmallEntries(C_hypre, 1e-15, 1);
+      C->WrapHypreParCSRMatrix(C_hypre);
+   }
+
+   if (job == BlockInverseScaleJob::RHS_ONLY ||
+       job == BlockInverseScaleJob::MATRIX_AND_RHS)
+   {
+      HypreParVector b_Hypre(A->GetComm(),
+                             A->GetGlobalNumRows(),
+                             b->GetData(), A->GetRowStarts());
+      hypre_ParVector *d_hypre;
+      hypre_ParvecBdiagInvScal(b_Hypre, blocksize, &d_hypre, *A);
+
+      d->WrapHypreParVector(d_hypre, true);
+   }
+}
+
+#endif
+
 #if MFEM_HYPRE_VERSION < 21400
 
 HypreParMatrix *Add(double alpha, const HypreParMatrix &A,
@@ -2817,6 +2905,11 @@ void HypreGMRES::SetTol(double tol)
    HYPRE_GMRESSetTol(gmres_solver, tol);
 }
 
+void HypreGMRES::SetAbsTol(double tol)
+{
+   HYPRE_GMRESSetAbsoluteTol(gmres_solver, tol);
+}
+
 void HypreGMRES::SetMaxIter(int max_iter)
 {
    HYPRE_GMRESSetMaxIter(gmres_solver, max_iter);
@@ -3620,6 +3713,116 @@ void HypreBoomerAMG::SetElasticityOptions(ParFiniteElementSpace *fespace)
    error_mode = IGNORE_HYPRE_ERRORS;
 }
 
+#if MFEM_HYPRE_VERSION >= 21800
+
+void HypreBoomerAMG::SetAdvectiveOptions(int distanceR,
+                                         const std::string &prerelax,
+                                         const std::string &postrelax)
+{
+   // Hypre parameters
+   int Sabs = 0;
+   int interp_type = 100;
+   int relax_type = 10;
+   int coarsen_type = 6;
+   double strength_tolC = 0.1;
+   double strength_tolR = 0.01;
+   double filter_tolR = 0.0;
+   double filterA_tol = 0.0;
+
+   // Set relaxation on specified grid points
+   int ns_down, ns_up, ns_coarse;
+   if (distanceR > 0)
+   {
+      ns_down = prerelax.length();
+      ns_up = postrelax.length();
+      ns_coarse = 1;
+
+      // Array to store relaxation scheme and pass to Hypre
+      HYPRE_Int **grid_relax_points = mfem_hypre_TAlloc(HYPRE_Int*, 4);
+      grid_relax_points[0] = NULL;
+      grid_relax_points[1] = mfem_hypre_TAlloc(HYPRE_Int, ns_down);
+      grid_relax_points[2] = mfem_hypre_TAlloc(HYPRE_Int, ns_up);
+      grid_relax_points[3] = mfem_hypre_TAlloc(HYPRE_Int, 1);
+      grid_relax_points[3][0] = 0;
+
+      // set down relax scheme
+      for (int i = 0; i<ns_down; i++)
+      {
+         if (prerelax[i] == 'F')
+         {
+            grid_relax_points[1][i] = -1;
+         }
+         else if (prerelax[i] == 'C')
+         {
+            grid_relax_points[1][i] = 1;
+         }
+         else if (prerelax[i] == 'A')
+         {
+            grid_relax_points[1][i] = 0;
+         }
+      }
+
+      // set up relax scheme
+      for (int i = 0; i<ns_up; i++)
+      {
+         if (postrelax[i] == 'F')
+         {
+            grid_relax_points[2][i] = -1;
+         }
+         else if (postrelax[i] == 'C')
+         {
+            grid_relax_points[2][i] = 1;
+         }
+         else if (postrelax[i] == 'A')
+         {
+            grid_relax_points[2][i] = 0;
+         }
+      }
+
+      HYPRE_BoomerAMGSetRestriction(amg_precond, distanceR);
+
+      HYPRE_BoomerAMGSetGridRelaxPoints(amg_precond, grid_relax_points);
+
+      HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type);
+   }
+
+   if (Sabs)
+   {
+      HYPRE_BoomerAMGSetSabs(amg_precond, Sabs);
+   }
+
+   HYPRE_BoomerAMGSetCoarsenType(amg_precond, coarsen_type);
+
+   // does not support aggressive coarsening
+   HYPRE_BoomerAMGSetAggNumLevels(amg_precond, 0);
+
+   HYPRE_BoomerAMGSetStrongThreshold(amg_precond, strength_tolC);
+
+   if (distanceR > 0)
+   {
+      HYPRE_BoomerAMGSetStrongThresholdR(amg_precond, strength_tolR);
+      HYPRE_BoomerAMGSetFilterThresholdR(amg_precond, filter_tolR);
+   }
+
+   if (relax_type > -1)
+   {
+      HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type);
+   }
+
+   if (distanceR > 0)
+   {
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_coarse, 3);
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_down,   1);
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_up,     2);
+
+      HYPRE_BoomerAMGSetADropTol(amg_precond, filterA_tol);
+      // type = -1: drop based on row inf-norm
+      HYPRE_BoomerAMGSetADropType(amg_precond, -1);
+   }
+}
+
+#endif
+
 HypreBoomerAMG::~HypreBoomerAMG()
 {
    for (int i = 0; i < rbms.Size(); i++)

linalg/hypre.hpp

@@ -81,6 +81,14 @@ private:
    inline void _SetDataAndSize_();
 
 public:
+
+   /// Default constructor, no underlying @a hypre_ParVector is created.
+   HypreParVector()
+   {
+      own_ParVector = false;
+      x = NULL;
+   }
+
    /** @brief Creates vector with given global size and parallel partitioning of
        the rows/columns given by @a col. */
    /** @anchor hypre_partitioning_descr
@@ -113,6 +121,9 @@ public:
    /// MPI communicator
    MPI_Comm GetComm() { return x->comm; }
 
+   /// Converts hypre's format to HypreParVector
+   void WrapHypreParVector(hypre_ParVector *y, bool owner=true);
+
    /// Returns the parallel row/column partitioning
    /** See @ref hypre_partitioning_descr "here" for a description of the
        partitioning array. */
@@ -228,15 +239,24 @@ public:
    /// An empty matrix to be used as a reference to an existing matrix
    HypreParMatrix();
 
+   /// Converts hypre's format to HypreParMatrix
+   /** If @a owner is false, ownership of @a a is not transferred */
+   void WrapHypreParCSRMatrix(hypre_ParCSRMatrix *a, bool owner = true)
+   {
+      Destroy();
+      Init();
+      A = a;
+      ParCSROwner = owner;
+      height = GetNumRows();
+      width = GetNumCols();
+   }
+
    /// Converts hypre's format to HypreParMatrix
    /** If @a owner is false, ownership of @a a is not transferred */
    explicit HypreParMatrix(hypre_ParCSRMatrix *a, bool owner = true)
    {
       Init();
-      A = a;
-      if (!owner) { ParCSROwner = 0; }
-      height = GetNumRows();
-      width = GetNumCols();
+      WrapHypreParCSRMatrix(a, owner);
    }
 
    /// Creates block-diagonal square parallel matrix.
@@ -387,6 +407,13 @@ public:
    void GetDiag(SparseMatrix &diag) const;
    /// Get the local off-diagonal block. NOTE: 'offd' will not own any data.
    void GetOffd(SparseMatrix &offd, HYPRE_Int* &cmap) const;
+   /** @brief Get a single SparseMatrix containing all rows from this processor,
+       merged from the diagonal and off-diagonal blocks stored by the
+       HypreParMatrix. */
+   /** @note The number of columns in the SparseMatrix will be the global number
+       of columns in the parallel matrix, so using this method may result in an
+       integer overflow in the column indices. */
+   void MergeDiagAndOffd(SparseMatrix &merged);
 
    /** Split the matrix into M x N equally sized blocks of parallel matrices.
        The size of 'blocks' must already be set to M x N. */
@@ -397,6 +424,13 @@ public:
    /// Returns the transpose of *this
    HypreParMatrix * Transpose() const;
 
+   /** Returns principle submatrix given by array of indices of connections
+       with relative size > @a threshold in *this. */
+#if MFEM_HYPRE_VERSION >= 21800
+   HypreParMatrix *ExtractSubmatrix(const Array<int> &indices,
+                                    double threshhold=0.0) const;
+#endif
+
    /// Returns the number of rows in the diagonal block of the ParCSRMatrix
    int GetNumRows() const
    {
@@ -549,6 +583,23 @@ public:
    Type GetType() const { return Hypre_ParCSR; }
 };
 
+#if MFEM_HYPRE_VERSION >= 21800
+
+enum class BlockInverseScaleJob
+{
+   MATRIX_ONLY,
+   RHS_ONLY,
+   MATRIX_AND_RHS
+};
+
+/** Constructs and applies block diagonal inverse of HypreParMatrix.
+    The enum @a job specifies whether the matrix or the RHS should be
+    scaled (or both). */
+void BlockInverseScale(const HypreParMatrix *A, HypreParMatrix *C,
+                       const Vector *b, HypreParVector *d,
+                       int blocksize, BlockInverseScaleJob job);
+#endif
+
 /** @brief Return a new matrix `C = alpha*A + beta*B`, assuming that both `A`
     and `B` use the same row and column partitions and the same `col_map_offd`
     arrays. */
@@ -649,11 +700,12 @@ public:
        4    = truncated l1-scaled block Gauss-Seidel/SSOR
        5    = lumped Jacobi
        6    = Gauss-Seidel
+       10   = On-processor forward solve for matrix w/ triangular structure
        16   = Chebyshev
        1001 = Taubin polynomial smoother
        1002 = FIR polynomial smoother. */
    enum Type { Jacobi = 0, l1Jacobi = 1, l1GS = 2, l1GStr = 4, lumpedJacobi = 5,
-               GS = 6, Chebyshev = 16, Taubin = 1001, FIR = 1002
+               GS = 6, OPFS = 10, Chebyshev = 16, Taubin = 1001, FIR = 1002
              };
 
    HypreSmoother();
@@ -767,6 +819,28 @@ public:
    virtual ~HypreSolver();
 };
 
+
+#if MFEM_HYPRE_VERSION >= 21800
+/** Preconditioner for HypreParMatrices that are triangular in some ordering.
+   Finds correct ordering and performs forward substitution on processor
+   as approximate inverse. Exact on one processor. */
+class HypreTriSolve : public HypreSolver
+{
+public:
+   HypreTriSolve() : HypreSolver() { }
+   explicit HypreTriSolve(HypreParMatrix &A) : HypreSolver(&A) { }
+   virtual operator HYPRE_Solver() const { return NULL; }
+
+   virtual HYPRE_PtrToParSolverFcn SetupFcn() const
+   { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSROnProcTriSetup; }
+   virtual HYPRE_PtrToParSolverFcn SolveFcn() const
+   { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSROnProcTriSolve; }
+
+   HypreParMatrix* GetData() { return A; }
+   virtual ~HypreTriSolve() { }
+};
+#endif
+
 /// PCG solver in hypre
 class HyprePCG : public HypreSolver
 {
@@ -844,6 +918,7 @@ public:
    virtual void SetOperator(const Operator &op);
 
    void SetTol(double tol);
+   void SetAbsTol(double tol);
    void SetMaxIter(int max_iter);
    void SetKDim(int dim);
    void SetLogging(int logging);
@@ -1128,9 +1203,94 @@ public:
        construct A. */
    void SetElasticityOptions(ParFiniteElementSpace *fespace);
 
+#if MFEM_HYPRE_VERSION >= 21800
+   /** Hypre parameters to use AIR AMG solve for advection-dominated problems.
+       See "Nonsymmetric Algebraic Multigrid Based on Local Approximate Ideal
+       Restriction (AIR)," Manteuffel, Ruge, Southworth, SISC (2018),
+       DOI:/10.1137/17M1144350. Options: "distanceR" -> distance of neighbor
+       DOFs to buld restriction operator; options include 1, 2, and 15 (1.5).
+       Strings "prerelax" and "postrelax" indicate points to relax on:
+       F = F-points, C = C-points, A = all points. E.g., FFC -> relax on
+       F-points, relax again on F-points, then relax on C-points. */
+   void SetAdvectiveOptions(int distance=15,  const std::string &prerelax="",
+                            const std::string &postrelax="FFC");
+
+   /// Expert option - consult hypre documentation/team
+   void SetStrongThresholdR(double strengthR)
+   { HYPRE_BoomerAMGSetStrongThresholdR(amg_precond, strengthR); }
+
+   /// Expert option - consult hypre documentation/team
+   void SetFilterThresholdR(double filterR)
+   { HYPRE_BoomerAMGSetFilterThresholdR(amg_precond, filterR); }
+
+   /// Expert option - consult hypre documentation/team
+   void SetRestriction(int restrict_type)
+   { HYPRE_BoomerAMGSetRestriction(amg_precond, restrict_type); }
+
+   /// Expert option - consult hypre documentation/team
+   void SetIsTriangular()
+   { HYPRE_BoomerAMGSetIsTriangular(amg_precond, 1); }
+
+   /// Expert option - consult hypre documentation/team
+   void SetGMRESSwitchR(int gmres_switch)
+   { HYPRE_BoomerAMGSetGMRESSwitchR(amg_precond, gmres_switch); }
+
+   /// Expert option - consult hypre documentation/team
+   void SetCycleNumSweeps(int prerelax, int postrelax)
+   {
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, prerelax,  1);
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, postrelax, 2);
+   }
+#endif
+
    void SetPrintLevel(int print_level)
    { HYPRE_BoomerAMGSetPrintLevel(amg_precond, print_level); }
 
+   void SetMaxIter(int max_iter)
+   { HYPRE_BoomerAMGSetMaxIter(amg_precond, max_iter); }
+
+   /// Expert option - consult hypre documentation/team
+   void SetMaxLevels(int max_levels)
+   { HYPRE_BoomerAMGSetMaxLevels(amg_precond, max_levels); }
+
+   /// Expert option - consult hypre documentation/team
+   void SetTol(double tol)
+   { HYPRE_BoomerAMGSetTol(amg_precond, tol); }
+
+   /// Expert option - consult hypre documentation/team
+   void SetStrengthThresh(double strength)
+   { HYPRE_BoomerAMGSetStrongThreshold(amg_precond, strength); }
+
+   /// Expert option - consult hypre documentation/team
+   void SetInterpolation(int interp_type)
+   { HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type); }
+
+   /// Expert option - consult hypre documentation/team
+   void SetCoarsening(int coarsen_type)
+   { HYPRE_BoomerAMGSetCoarsenType(amg_precond, coarsen_type); }
+
+   /// Expert option - consult hypre documentation/team
+   void SetRelaxType(int relax_type)
+   { HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type); }
+
+   /// Expert option - consult hypre documentation/team
+   void SetCycleType(int cycle_type)
+   { HYPRE_BoomerAMGSetCycleType(amg_precond, cycle_type); }
+
+   void GetNumIterations(int &num_it)
+   { HYPRE_BoomerAMGGetNumIterations(amg_precond, &num_it); }
+
+   /// Expert option - consult hypre documentation/team
+   void SetNodal(int blocksize)
+   {
+      HYPRE_BoomerAMGSetNumFunctions(amg_precond, blocksize);
+      HYPRE_BoomerAMGSetNodal(amg_precond, 1);
+   }
+
+   /// Expert option - consult hypre documentation/team
+   void SetAggressiveCoarsening(int num_levels)
+   { HYPRE_BoomerAMGSetAggNumLevels(amg_precond, num_levels); }
+
    /// The typecast to HYPRE_Solver returns the internal amg_precond
    virtual operator HYPRE_Solver() const { return amg_precond; }
 
@@ -1139,6 +1299,8 @@ public:
    virtual HYPRE_PtrToParSolverFcn SolveFcn() const
    { return (HYPRE_PtrToParSolverFcn) HYPRE_BoomerAMGSolve; }
 
+   using HypreSolver::Mult;
+
    virtual ~HypreBoomerAMG();
 };
 

makefile

@@ -684,9 +684,10 @@ status info:
 ASTYLE_BIN = astyle
 ASTYLE = $(ASTYLE_BIN) --options=$(SRC)config/mfem.astylerc
 ASTYLE_VER = "Artistic Style Version 2.05.1"
-FORMAT_FILES = $(foreach dir,$(DIRS) $(EM_DIRS) config,"$(dir)/*.?pp")
-FORMAT_FILES += "tests/unit/*.cpp"
-FORMAT_FILES += $(foreach dir,general linalg mesh fem,"tests/unit/$(dir)/*.?pp")
+FORMAT_FILES := $(foreach dir,$(DIRS) $(EM_DIRS) config,$(dir)/*.?pp)
+FORMAT_FILES += tests/unit/*.cpp
+FORMAT_FILES += $(foreach dir,general linalg mesh fem,tests/unit/$(dir)/*.?pp)
+FORMAT_FILES := $(filter-out general/tinyxml2.cpp,$(wildcard $(FORMAT_FILES)))
 
 COUT_CERR_FILES = $(foreach dir,$(DIRS),$(dir)/*.[ch]pp)
 COUT_CERR_EXCLUDE = '^general/error\.cpp' '^general/globals\.[ch]pp'

mesh/mesh.cpp

@@ -3448,6 +3448,7 @@ void Mesh::Loader(std::istream &input, int generate_edges,
 
    Clear();
 
+   istream::pos_type beginning_pos = input.tellg();
    string mesh_type;
    input >> ws;
    getline(input, mesh_type);
@@ -3489,11 +3490,20 @@ void Mesh::Loader(std::istream &input, int generate_edges,
    {
       ReadTrueGridMesh(input);
    }
-   else if (mesh_type == "# vtk DataFile Version 3.0" ||
-            mesh_type == "# vtk DataFile Version 2.0") // VTK
+   else if (mesh_type.rfind("# vtk DataFile Version") == 0)
    {
+      int major_vtk_version = mesh_type[mesh_type.length()-3] - '0';
+      // int minor_vtk_version = mesh_type[mesh_type.length()-1] - '0';
+      MFEM_VERIFY(major_vtk_version >= 2 && major_vtk_version <= 4,
+                  "Unsupported VTK format");
       ReadVTKMesh(input, curved, read_gf, finalize_topo);
    }
+   else if (mesh_type.rfind("<VTKFile ") == 0)
+   {
+      // Go back to beginning of stream
+      input.seekg(beginning_pos);
+      ReadXML_VTKMesh(input, curved, read_gf, finalize_topo);
+   }
    else if (mesh_type == "MFEM NURBS mesh v1.0")
    {
       ReadNURBSMesh(input, curved, read_gf);
@@ -8938,7 +8948,7 @@ void Mesh::PrintVTK(std::ostream &out)
          const int nv = elements[i]->GetNVertices();
          out << nv;
          Geometry::Type geom = elements[i]->GetGeometryType();
-         const int *perm = (geom == Geometry::PRISM) ? vtk_prism_perm : NULL;
+         const int *perm = VTKGeometry::VertexPermutation[geom];
          for (int j = 0; j < nv; j++)
          {
             out << ' ' << v[perm ? perm[j] : j];
@@ -9023,35 +9033,9 @@ void Mesh::PrintVTK(std::ostream &out)
    for (int i = 0; i < NumOfElements; i++)
    {
       int vtk_cell_type = 5;
-      Geometry::Type geom_type = GetElement(i)->GetGeometryType();
-      if (order == 1)
-      {
-         switch (geom_type)
-         {
-            case Geometry::POINT:        vtk_cell_type = 1;   break;
-            case Geometry::SEGMENT:      vtk_cell_type = 3;   break;
-            case Geometry::TRIANGLE:     vtk_cell_type = 5;   break;
-            case Geometry::SQUARE:       vtk_cell_type = 9;   break;
-            case Geometry::TETRAHEDRON:  vtk_cell_type = 10;  break;
-            case Geometry::CUBE:         vtk_cell_type = 12;  break;
-            case Geometry::PRISM:        vtk_cell_type = 13;  break;
-            default: break;
-         }
-      }
-      else if (order == 2)
-      {
-         switch (geom_type)
-         {
-            case Geometry::SEGMENT:      vtk_cell_type = 21;  break;
-            case Geometry::TRIANGLE:     vtk_cell_type = 22;  break;
-            case Geometry::SQUARE:       vtk_cell_type = 28;  break;
-            case Geometry::TETRAHEDRON:  vtk_cell_type = 24;  break;
-            case Geometry::CUBE:         vtk_cell_type = 29;  break;
-            case Geometry::PRISM:        vtk_cell_type = 32;  break;
-            default: break;
-         }
-      }
-
+      Geometry::Type geom = GetElement(i)->GetGeometryType();
+      if (order == 1) { vtk_cell_type = VTKGeometry::Map[geom]; }
+      else if (order == 2) { vtk_cell_type = VTKGeometry::QuadraticMap[geom]; }
       out << vtk_cell_type << '\n';
    }
 
@@ -9269,7 +9253,7 @@ void Mesh::PrintVTU(std::ostream &out, int ref, VTKFormat format,
          {
             coff = coff+nv;
             offset.push_back(coff);
-            const int *p = (geom == Geometry::PRISM) ? vtk_prism_perm : NULL;
+            const int *p = VTKGeometry::VertexPermutation[geom];
             for (int k = 0; k < nv; k++, j++)
             {
                WriteBinaryOrASCII(out, buf, np + RG[p ? p[j] : j], " ", format);
@@ -9300,39 +9284,14 @@ void Mesh::PrintVTU(std::ostream &out, int ref, VTKFormat format,
    out << "<DataArray type=\"UInt8\" Name=\"types\" format=\""
        << fmt_str << "\">" << std::endl;
    // cell types
+   const int *vtk_geom_map =
+      high_order_output ? VTKGeometry::HighOrderMap : VTKGeometry::Map;
    for (int i = 0; i < ne; i++)
    {
       Geometry::Type geom = get_geom(i);
       uint8_t vtk_cell_type = 5;
 
-      // VTK element types defined at: https://git.io/JvZLm
-      switch (geom)
-      {
-         case Geometry::POINT:
-            vtk_cell_type = 1;
-            break;
-         case Geometry::SEGMENT:
-            vtk_cell_type = high_order_output ? 68 : 3;
-            break;
-         case Geometry::TRIANGLE:
-            vtk_cell_type = high_order_output ? 69 : 5;
-            break;
-         case Geometry::SQUARE:
-            vtk_cell_type = high_order_output ? 70 : 9;
-            break;
-         case Geometry::TETRAHEDRON:
-            vtk_cell_type = high_order_output ? 71 : 10;
-            break;
-         case Geometry::CUBE:
-            vtk_cell_type = high_order_output ? 72 : 12;
-            break;
-         case Geometry::PRISM:
-            vtk_cell_type = high_order_output ? 73 : 13;
-            break;
-         default:
-            MFEM_ABORT("Unrecognized VTK element type \"" << geom << "\"");
-            break;
-      }
+      vtk_cell_type = vtk_geom_map[geom];
 
       if (high_order_output)
       {
@@ -9481,21 +9440,7 @@ void Mesh::PrintVTK(std::ostream &out, int ref, int field_data)
       int nv = Geometries.GetVertices(geom)->GetNPoints();
       RefG = GlobGeometryRefiner.Refine(geom, ref, 1);
       Array<int> &RG = RefG->RefGeoms;
-      int vtk_cell_type = 5;
-
-      switch (geom)
-      {
-         case Geometry::POINT:        vtk_cell_type = 1;   break;
-         case Geometry::SEGMENT:      vtk_cell_type = 3;   break;
-         case Geometry::TRIANGLE:     vtk_cell_type = 5;   break;
-         case Geometry::SQUARE:       vtk_cell_type = 9;   break;
-         case Geometry::TETRAHEDRON:  vtk_cell_type = 10;  break;
-         case Geometry::CUBE:         vtk_cell_type = 12;  break;
-         case Geometry::PRISM:        vtk_cell_type = 13;  break;
-         default:
-            MFEM_ABORT("Unrecognized VTK element type \"" << geom << "\"");
-            break;
-      }
+      int vtk_cell_type = VTKGeometry::Map[geom];
 
       for (int j = 0; j < RG.Size(); j += nv)
       {

mesh/mesh.hpp

@@ -239,8 +239,15 @@ protected:
    void ReadNetgen2DMesh(std::istream &input, int &curved);
    void ReadNetgen3DMesh(std::istream &input);
    void ReadTrueGridMesh(std::istream &input);
+   void CreateVTKMesh(const Vector &points, const Array<int> &cell_data,
+                      const Array<int> &cell_offsets,
+                      const Array<int> &cell_types,
+                      const Array<int> &cell_attributes,
+                      int &curved, int &read_gf, bool &finalize_topo);
    void ReadVTKMesh(std::istream &input, int &curved, int &read_gf,
                     bool &finalize_topo);
+   void ReadXML_VTKMesh(std::istream &input, int &curved, int &read_gf,
+                        bool &finalize_topo);
    void ReadNURBSMesh(std::istream &input, int &curved, int &read_gf);
    void ReadInlineMesh(std::istream &input, bool generate_edges = false);
    void ReadGmshMesh(std::istream &input, int &curved, int &read_gf);

mesh/mesh_readers.cpp

@@ -12,10 +12,13 @@
 #include "mesh_headers.hpp"
 #include "../fem/fem.hpp"
 #include "../general/text.hpp"
+#include "../general/tinyxml2.h"
 #include "gmsh.hpp"
 
 #include <iostream>
 #include <cstdio>
+#include <vector>
+#include <algorithm>
 
 #ifdef MFEM_USE_NETCDF
 #include "netcdf.h"
@@ -373,6 +376,398 @@ const int Mesh::vtk_quadratic_hex[27] =
    24, 22, 21, 23, 20, 25, 26
 };
 
+void Mesh::CreateVTKMesh(const Vector &points, const Array<int> &cell_data,
+                         const Array<int> &cell_offsets,
+                         const Array<int> &cell_types,
+                         const Array<int> &cell_attributes,
+                         int &curved, int &read_gf, bool &finalize_topo)
+{
+   int np = points.Size()/3;
+   Dim = -1;
+   NumOfElements = cell_types.Size();
+   elements.SetSize(NumOfElements);
+
+   int order = -1;
+   bool legacy_elem = false, lagrange_elem = false;
+
+   int j = 0;
+   for (int i = 0; i < NumOfElements; i++)
+   {
+      int ct = cell_types[i];
+      Geometry::Type geom = VTKGeometry::GetMFEMGeometry(ct);
+      elements[i] = NewElement(geom);
+      if (cell_attributes.Size() > 0)
+      {
+         elements[i]->SetAttribute(cell_attributes[i]);
+      }
+      // VTK ordering of vertices is the same as MFEM ordering of vertices
+      // for all element types *except* prisms, which require a permutation
+      if (geom == Geometry::PRISM && ct != VTKGeometry::LAGRANGE_PRISM)
+      {
+         int prism_vertices[6];
+         for (int k=0; k<6; ++k)
+         {
+            prism_vertices[k] = cell_data[j+VTKGeometry::PrismMap[k]];
+         }
+         elements[i]->SetVertices(prism_vertices);
+      }
+      else
+      {
+         elements[i]->SetVertices(&cell_data[j]);
+      }
+
+      int elem_dim = Geometry::Dimension[geom];
+      int elem_order = VTKGeometry::GetOrder(ct, cell_offsets[i] - j);
+
+      if (VTKGeometry::IsLagrange(ct)) { lagrange_elem = true; }
+      else { legacy_elem = true; }
+
+      MFEM_VERIFY(Dim == -1 || Dim == elem_dim,
+                  "Elements with different dimensions are not supported");
+      MFEM_VERIFY(order == -1 || order == elem_order,
+                  "Elements with different orders are not supported");
+      MFEM_VERIFY(legacy_elem != lagrange_elem,
+                  "Mixing of legacy and Lagrange cell types is not supported");
+      Dim = elem_dim;
+      order = elem_order;
+      j = cell_offsets[i];
+   }
+
+   if (order == 1 && !lagrange_elem)
+   {
+      NumOfVertices = np;
+      vertices.SetSize(np);
+      for (int i = 0; i < np; i++)
+      {
+         vertices[i](0) = points(3*i+0);
+         vertices[i](1) = points(3*i+1);
+         vertices[i](2) = points(3*i+2);
+      }
+      // No boundary is defined in a VTK mesh
+      NumOfBdrElements = 0;
+      FinalizeTopology();
+      CheckElementOrientation(true);
+   }
+   else
+   {
+      // The following section of code is shared for legacy quadratic and the
+      // Lagrange high order elements
+      curved = 1;
+
+      // generate new enumeration for the vertices
+      Array<int> pts_dof(np);
+      pts_dof = -1;
+      // mark vertex points
+      for (int i = 0; i < NumOfElements; i++)
+      {
+         int *v = elements[i]->GetVertices();
+         int nv = elements[i]->GetNVertices();
+         for (int j = 0; j < nv; j++)
+         {
+            if (pts_dof[v[j]] == -1) { pts_dof[v[j]] = 0; }
+         }
+      }
+
+      // The following loop reorders pts_dofs so vertices are visited in
+      // canonical order
+
+      // Keep the original ordering of the vertices
+      int i, n;
+      for (n = i = 0; i < np; i++)
+      {
+         if (pts_dof[i] != -1)
+         {
+            pts_dof[i] = n++;
+         }
+      }
+      // update the element vertices
+      for (int i = 0; i < NumOfElements; i++)
+      {
+         int *v = elements[i]->GetVertices();
+         int nv = elements[i]->GetNVertices();
+         for (int j = 0; j < nv; j++)
+         {
+            v[j] = pts_dof[v[j]];
+         }
+      }
+      // Define the 'vertices' from the 'points' through the 'pts_dof' map
+      NumOfVertices = n;
+      vertices.SetSize(n);
+      for (int i = 0; i < np; i++)
+      {
+         int j = pts_dof[i];
+         if (j != -1)
+         {
+            vertices[j](0) = points(3*i+0);
+            vertices[j](1) = points(3*i+1);
+            vertices[j](2) = points(3*i+2);
+         }
+      }
+
+      // No boundary is defined in a VTK mesh
+      NumOfBdrElements = 0;
+
+      // determine spaceDim based on min/max differences detected each dimension
+      if (vertices.Size() > 0)
+      {
+         double min_value, max_value;
+         for (int d=0; d<3; ++d)
+         {
+            min_value = max_value = vertices[0](d);
+            for (int i = 1; i < vertices.Size(); i++)
+            {
+               min_value = std::min(min_value,vertices[i](d));
+               max_value = std::max(max_value,vertices[i](d));
+               if (min_value != max_value)
+               {
+                  spaceDim++;
+                  break;
+               }
+            }
+         }
+      }
+      // Generate faces and edges so that we can define
+      // FE space on the mesh
+      FinalizeTopology();
+
+      FiniteElementCollection *fec;
+      FiniteElementSpace *fes;
+      if (legacy_elem)
+      {
+         // Define quadratic FE space
+         fec = new QuadraticFECollection;
+         fes = new FiniteElementSpace(this, fec, spaceDim);
+         Nodes = new GridFunction(fes);
+         Nodes->MakeOwner(fec); // Nodes will destroy 'fec' and 'fes'
+         own_nodes = 1;
+
+         // Map vtk points to edge/face/element dofs
+         Array<int> dofs;
+         for (int i = 0; i < NumOfElements; i++)
+         {
+            fes->GetElementDofs(i, dofs);
+            const int *vtk_mfem;
+            switch (elements[i]->GetGeometryType())
+            {
+               case Geometry::TRIANGLE:
+               case Geometry::SQUARE:
+                  vtk_mfem = vtk_quadratic_hex; break; // identity map
+               case Geometry::TETRAHEDRON:
+                  vtk_mfem = vtk_quadratic_tet; break;
+               case Geometry::CUBE:
+                  vtk_mfem = vtk_quadratic_hex; break;
+               case Geometry::PRISM:
+                  vtk_mfem = vtk_quadratic_wedge; break;
+               default:
+                  vtk_mfem = NULL; // suppress a warning
+                  break;
+            }
+
+            int offset = (i == 0) ? 0 : cell_offsets[i-1];
+            for (int j = 0; j < dofs.Size(); j++)
+            {
+               if (pts_dof[cell_data[offset+j]] == -1)
+               {
+                  pts_dof[cell_data[offset+j]] = dofs[vtk_mfem[j]];
+               }
+               else
+               {
+                  if (pts_dof[cell_data[offset+j]] != dofs[vtk_mfem[j]])
+                  {
+                     MFEM_ABORT("VTK mesh: inconsistent quadratic mesh!");
+                  }
+               }
+            }
+         }
+      }
+      else
+      {
+         // Define H1 FE space
+         fec = new H1_FECollection(order,Dim,BasisType::ClosedUniform);
+         fes = new FiniteElementSpace(this, fec, spaceDim);
+         Nodes = new GridFunction(fes);
+         Nodes->MakeOwner(fec); // Nodes will destroy 'fec' and 'fes'
+         own_nodes = 1;
+         Array<int> dofs;
+
+         std::map<Geometry::Type,Array<int>> vtk_inv_maps;
+         std::map<Geometry::Type,const Array<int>*> lex_orderings;
+
+         int i, n;
+         for (n = i = 0; i < NumOfElements; i++)
+         {
+            Geometry::Type geom = GetElementBaseGeometry(i);
+            fes->GetElementDofs(i, dofs);
+
+            Array<int> &vtk_inv_map = vtk_inv_maps[geom];
+            if (vtk_inv_map.Size() == 0)
+            {
+               Array<int> vtk_map;
+               CreateVTKElementConnectivity(vtk_map, geom, order);
+               vtk_inv_map.SetSize(vtk_map.Size());
+               for (int j=0; j<vtk_map.Size(); ++j)
+               {
+                  vtk_inv_map[vtk_map[j]] = j;
+               }
+            }
+            const Array<int> *&lex_ordering = lex_orderings[geom];
+            if (!lex_ordering)
+            {
+               const FiniteElement *fe = fes->GetFE(i);
+               const NodalFiniteElement *nodal_fe =
+                  dynamic_cast<const NodalFiniteElement*>(fe);
+               MFEM_ASSERT(nodal_fe != NULL, "Unsupported element type");
+               lex_ordering = &nodal_fe->GetLexicographicOrdering();
+            }
+
+            for (int lex_idx = 0; lex_idx < dofs.Size(); lex_idx++)
+            {
+               int mfem_idx = (*lex_ordering)[lex_idx];
+               int vtk_idx = vtk_inv_map[lex_idx];
+               int pt_idx = cell_data[n + vtk_idx];
+               if (pts_dof[pt_idx] == -1)
+               {
+                  pts_dof[pt_idx] = dofs[mfem_idx];
+               }
+               else
+               {
+                  if (pts_dof[pt_idx] != dofs[mfem_idx])
+                  {
+                     MFEM_ABORT("VTK mesh: inconsistent Lagrange mesh!");
+                  }
+               }
+            }
+            n += dofs.Size();
+         }
+      }
+      // Define the 'Nodes' from the 'points' through the 'pts_dof' map
+      Array<int> dofs;
+      for (int i = 0; i < np; i++)
+      {
+         dofs.SetSize(1);
+         if (pts_dof[i] != -1)
+         {
+            dofs[0] = pts_dof[i];
+            fes->DofsToVDofs(dofs);
+            for (int d = 0; d < dofs.Size(); d++)
+            {
+               (*Nodes)(dofs[d]) = points(3*i+d);
+            }
+         }
+      }
+      read_gf = 0;
+   }
+}
+
+void Mesh::ReadXML_VTKMesh(std::istream &input, int &curved, int &read_gf,
+                           bool &finalize_topo)
+{
+   using namespace tinyxml2;
+
+   const char *erstr = "XML parsing error";
+
+   // Read entire stream into buffer
+   std::istreambuf_iterator<char> eos;
+   std::vector<char> buf(std::istreambuf_iterator<char>(input), eos);
+   buf.push_back('\0'); // null-terminate buffer
+
+   XMLDocument xml;
+   xml.Parse(buf.data());
+   MFEM_VERIFY(xml.ErrorID() == XML_SUCCESS, erstr);
+
+   const XMLElement *vtkfile = xml.FirstChildElement();
+   MFEM_VERIFY(vtkfile, erstr);
+   MFEM_VERIFY(std::string(vtkfile->Name()) == "VTKFile", erstr);
+   const XMLElement *vtu = vtkfile->FirstChildElement();
+   MFEM_VERIFY(vtu, erstr);
+   MFEM_VERIFY(std::string(vtu->Name()) == "UnstructuredGrid", erstr);
+
+   // Count the number of points and cells
+   const XMLElement *piece = vtu->FirstChildElement();
+   MFEM_VERIFY(std::string(piece->Name()) == "Piece", erstr);
+   MFEM_VERIFY(piece->NextSiblingElement() == NULL,
+               "XML VTK meshes with more than one Piece are not supported");
+   int npts = piece->IntAttribute("NumberOfPoints");
+   int ncells = piece->IntAttribute("NumberOfCells");
+
+   // Read the points
+   Vector points(3*npts);
+   const XMLElement *pts_xml;
+   for (pts_xml = piece->FirstChildElement();
+        pts_xml != NULL;
+        pts_xml = pts_xml->NextSiblingElement())
+   {
+      if (std::string(pts_xml->Name()) == "Points")
+      {
+         const XMLElement *pts_data = pts_xml->FirstChildElement();
+         MFEM_VERIFY(std::string(pts_data->Name()) == "DataArray", erstr);
+         MFEM_VERIFY(std::string(pts_data->Attribute("Name")) == "Points",
+                     erstr);
+         MFEM_VERIFY(pts_data->IntAttribute("NumberOfComponents") == 3,
+                     "XML VTK Points DataArray must have 3 components");
+         const char *pts_txt = pts_data->GetText();
+         MFEM_VERIFY(pts_txt != NULL, erstr);
+         std::istringstream pts_stream(pts_txt);
+         points.Load(pts_stream, 3*npts);
+         break;
+      }
+   }
+   if (pts_xml == NULL) { MFEM_ABORT(erstr); }
+
+   // Read the cells
+   Array<int> cell_data, cell_offsets, cell_types;
+   const XMLElement *cells_xml;
+   for (cells_xml = piece->FirstChildElement();
+        cells_xml != NULL;
+        cells_xml = cells_xml->NextSiblingElement())
+   {
+      if (std::string(cells_xml->Name()) == "Cells")
+      {
+         const char *cell_data_txt = NULL;
+         for (const XMLElement *data_xml = cells_xml->FirstChildElement();
+              data_xml != NULL;
+              data_xml = data_xml->NextSiblingElement())
+         {
+            MFEM_VERIFY(std::string(data_xml->Name()) == "DataArray", erstr);
+            std::string data_name(data_xml->Attribute("Name"));
+            const char *data_txt = data_xml->GetText();
+            MFEM_VERIFY(data_txt != NULL, erstr);
+            if (data_name == "offsets")
+            {
+               std::istringstream data_stream(data_txt);
+               cell_offsets.Load(ncells, data_stream);
+            }
+            else if (data_name == "types")
+            {
+               std::istringstream data_stream(data_txt);
+               cell_types.Load(ncells, data_stream);
+            }
+            else if (data_name == "connectivity")
+            {
+               // Have to read the connectivity after the offsets, because we
+               // don't know how many points to read until we have the offsets
+               // (size of connectivity array is equal to the last offset), so
+               // store the data pointer and read this array later.
+               cell_data_txt = data_txt;
+            }
+         }
+         MFEM_VERIFY(cell_offsets.Size() == ncells, erstr);
+         MFEM_VERIFY(cell_types.Size() == ncells, erstr);
+         MFEM_VERIFY(cell_data_txt != NULL, erstr);
+         int cell_data_size = cell_offsets.Last();
+         std::istringstream cell_data_stream(cell_data_txt);
+         cell_data.Load(cell_data_size, cell_data_stream);
+         break;
+      }
+   }
+   if (cells_xml == NULL) { MFEM_ABORT(erstr); }
+
+   // Currently don't support reading cell attributes from VTK mesh
+   Array<int> cell_attributes;
+   CreateVTKMesh(points, cell_data, cell_offsets, cell_types, cell_attributes,
+                 curved, read_gf, finalize_topo);
+}
+
 void Mesh::ReadVTKMesh(std::istream &input, int &curved, int &read_gf,
                        bool &finalize_topo)
 {
@@ -382,8 +777,6 @@ void Mesh::ReadVTKMesh(std::istream &input, int &curved, int &read_gf,
    //   * https://lorensen.github.io/VTKExamples/site/VTKFileFormats
    //   * https://www.kitware.com/products/books/VTKUsersGuide.pdf
 
-   int i, j, n, attr;
-
    string buff;
    getline(input, buff); // comment line
    getline(input, buff);
@@ -412,141 +805,78 @@ void Mesh::ReadVTKMesh(std::istream &input, int &curved, int &read_gf,
       }
    }
    while (buff != "POINTS");
-   int np = 0;
+
    Vector points;
+   int np;
+   input >> np >> ws;
+   getline(input, buff); // "double"
+   points.Load(input, 3*np);
+
+   //skip metadata
+   // Looks like:
+   // METADATA
+   //INFORMATION 2
+   //NAME L2_NORM_RANGE LOCATION vtkDataArray
+   //DATA 2 0 5.19615
+   //NAME L2_NORM_FINITE_RANGE LOCATION vtkDataArray
+   //DATA 2 0 5.19615
+   do
    {
-      input >> np >> ws;
-      points.SetSize(3*np);
-      getline(input, buff); // "double"
-      for (i = 0; i < points.Size(); i++)
+      input >> buff;
+      if (!input.good())
       {
-         input >> points(i);
+         MFEM_ABORT("VTK mesh does not have CELLS data!");
       }
    }
+   while (buff != "CELLS");
 
    // Read the cells
-   NumOfElements = n = 0;
-   Array<int> cells_data;
-   input >> ws >> buff;
+   Array<int> cell_data, cell_offsets;
    if (buff == "CELLS")
    {
-      input >> NumOfElements >> n >> ws;
-      cells_data.SetSize(n);
-      for (i = 0; i < n; i++)
+      int ncells, n;
+      input >> ncells >> n >> ws;
+      cell_offsets.SetSize(ncells);
+      cell_data.SetSize(n - ncells);
+      int offset = 0;
+      for (int i=0; i<ncells; ++i)
       {
-         input >> cells_data[i];
+         int nv;
+         input >> nv;
+         cell_offsets[i] = offset + nv;
+         for (int j=0; j<nv; ++j)
+         {
+            input >> cell_data[offset + j];
+         }
+         offset += nv;
       }
    }
 
    // Read the cell types
-   Dim = -1;
-   int order = -1;
    input >> ws >> buff;
+   Array<int> cell_types;
+   int ncells;
    if (buff == "CELL_TYPES")
    {
-      input >> NumOfElements;
-      elements.SetSize(NumOfElements);
-      for (j = i = 0; i < NumOfElements; i++)
-      {
-         int ct, elem_dim, elem_order = 1;
-         input >> ct;
-         switch (ct)
-         {
-            case 5:   // triangle
-               elem_dim = 2;
-               elements[i] = new Triangle(&cells_data[j+1]);
-               break;
-            case 9:   // quadrilateral
-               elem_dim = 2;
-               elements[i] = new Quadrilateral(&cells_data[j+1]);
-               break;
-            case 10:  // tetrahedron
-               elem_dim = 3;
-#ifdef MFEM_USE_MEMALLOC
-               elements[i] = TetMemory.Alloc();
-               elements[i]->SetVertices(&cells_data[j+1]);
-#else
-               elements[i] = new Tetrahedron(&cells_data[j+1]);
-#endif
-               break;
-            case 12:  // hexahedron
-               elem_dim = 3;
-               elements[i] = new Hexahedron(&cells_data[j+1]);
-               break;
-            case 13:  // wedge
-               elem_dim = 3;
-               // switch between vtk vertex ordering and mfem vertex ordering:
-               // swap vertices (1,2) and (4,5)
-               elements[i] =
-                  new Wedge(cells_data[j+1], cells_data[j+3], cells_data[j+2],
-                            cells_data[j+4], cells_data[j+6], cells_data[j+5]);
-               break;
-
-            case 22:  // quadratic triangle
-               elem_dim = 2;
-               elem_order = 2;
-               elements[i] = new Triangle(&cells_data[j+1]);
-               break;
-            case 28:  // biquadratic quadrilateral
-               elem_dim = 2;
-               elem_order = 2;
-               elements[i] = new Quadrilateral(&cells_data[j+1]);
-               break;
-            case 24:  // quadratic tetrahedron
-               elem_dim = 3;
-               elem_order = 2;
-#ifdef MFEM_USE_MEMALLOC
-               elements[i] = TetMemory.Alloc();
-               elements[i]->SetVertices(&cells_data[j+1]);
-#else
-               elements[i] = new Tetrahedron(&cells_data[j+1]);
-#endif
-               break;
-            case 32: // biquadratic-quadratic wedge
-               elem_dim = 3;
-               elem_order = 2;
-               // switch between vtk vertex ordering and mfem vertex ordering:
-               // swap vertices (1,2) and (4,5)
-               elements[i] =
-                  new Wedge(cells_data[j+1], cells_data[j+3], cells_data[j+2],
-                            cells_data[j+4], cells_data[j+6], cells_data[j+5]);
-               break;
-            case 29:  // triquadratic hexahedron
-               elem_dim = 3;
-               elem_order = 2;
-               elements[i] = new Hexahedron(&cells_data[j+1]);
-               break;
-            default:
-               MFEM_ABORT("VTK mesh : cell type " << ct << " is not supported!");
-               return;
-         }
-         MFEM_VERIFY(Dim == -1 || Dim == elem_dim,
-                     "elements with different dimensions are not supported");
-         MFEM_VERIFY(order == -1 || order == elem_order,
-                     "elements with different orders are not supported");
-         Dim = elem_dim;
-         order = elem_order;
-         j += cells_data[j] + 1;
-      }
+      input >> ncells;
+      cell_types.Load(ncells, input);
    }
 
-   // Read attributes
+   // Read cell attributes
    streampos sp = input.tellg();
    input >> ws >> buff;
+   Array<int> cell_attributes;
    if (buff == "CELL_DATA")
    {
+      int n;
       input >> n >> ws;
       getline(input, buff);
       filter_dos(buff);
       // "SCALARS material dataType numComp"
-      if (!strncmp(buff.c_str(), "SCALARS material", 16))
+      if (buff.rfind("SCALARS material") == 0)
       {
          getline(input, buff); // "LOOKUP_TABLE default"
-         for (i = 0; i < NumOfElements; i++)
-         {
-            input >> attr;
-            elements[i]->SetAttribute(attr);
-         }
+         cell_attributes.Load(ncells, input);
       }
       else
       {
@@ -558,138 +888,9 @@ void Mesh::ReadVTKMesh(std::istream &input, int &curved, int &read_gf,
       input.seekg(sp);
    }
 
-   if (order == 1)
-   {
-      cells_data.DeleteAll();
-      NumOfVertices = np;
-      vertices.SetSize(np);
-      for (i = 0; i < np; i++)
-      {
-         vertices[i](0) = points(3*i+0);
-         vertices[i](1) = points(3*i+1);
-         vertices[i](2) = points(3*i+2);
-      }
-      points.Destroy();
-
-      // No boundary is defined in a VTK mesh
-      NumOfBdrElements = 0;
-   }
-   else if (order == 2)
-   {
-      curved = 1;
-
-      // generate new enumeration for the vertices
-      Array<int> pts_dof(np);
-      pts_dof = -1;
-      for (n = i = 0; i < NumOfElements; i++)
-      {
-         int *v = elements[i]->GetVertices();
-         int nv = elements[i]->GetNVertices();
-         for (j = 0; j < nv; j++)
-            if (pts_dof[v[j]] == -1)
-            {
-               pts_dof[v[j]] = n++;
-            }
-      }
-      // keep the original ordering of the vertices
-      for (n = i = 0; i < np; i++)
-         if (pts_dof[i] != -1)
-         {
-            pts_dof[i] = n++;
-         }
-      // update the element vertices
-      for (i = 0; i < NumOfElements; i++)
-      {
-         int *v = elements[i]->GetVertices();
-         int nv = elements[i]->GetNVertices();
-         for (j = 0; j < nv; j++)
-         {
-            v[j] = pts_dof[v[j]];
-         }
-      }
-      // Define the 'vertices' from the 'points' through the 'pts_dof' map
-      NumOfVertices = n;
-      vertices.SetSize(n);
-      for (i = 0; i < np; i++)
-      {
-         if ((j = pts_dof[i]) != -1)
-         {
-            vertices[j](0) = points(3*i+0);
-            vertices[j](1) = points(3*i+1);
-            vertices[j](2) = points(3*i+2);
-         }
-      }
-
-      // No boundary is defined in a VTK mesh
-      NumOfBdrElements = 0;
-
-      // Generate faces and edges so that we can define quadratic
-      // FE space on the mesh
-      FinalizeTopology();
-      finalize_topo = false;
-
-      // Define quadratic FE space
-      FiniteElementCollection *fec = new QuadraticFECollection;
-      FiniteElementSpace *fes = new FiniteElementSpace(this, fec, Dim);
-      Nodes = new GridFunction(fes);
-      Nodes->MakeOwner(fec); // Nodes will destroy 'fec' and 'fes'
-      own_nodes = 1;
-
-      // Map vtk points to edge/face/element dofs
-      Array<int> dofs;
-      for (n = i = 0; i < NumOfElements; i++)
-      {
-         fes->GetElementDofs(i, dofs);
-         const int *vtk_mfem;
-         switch (elements[i]->GetGeometryType())
-         {
-            case Geometry::TRIANGLE:
-            case Geometry::SQUARE:
-               vtk_mfem = vtk_quadratic_hex; break; // identity map
-            case Geometry::TETRAHEDRON:
-               vtk_mfem = vtk_quadratic_tet; break;
-            case Geometry::CUBE:
-               vtk_mfem = vtk_quadratic_hex; break;
-            case Geometry::PRISM:
-               vtk_mfem = vtk_quadratic_wedge; break;
-            default:
-               vtk_mfem = NULL; // suppress a warning
-               break;
-         }
-
-         for (n++, j = 0; j < dofs.Size(); j++, n++)
-         {
-            if (pts_dof[cells_data[n]] == -1)
-            {
-               pts_dof[cells_data[n]] = dofs[vtk_mfem[j]];
-            }
-            else
-            {
-               if (pts_dof[cells_data[n]] != dofs[vtk_mfem[j]])
-               {
-                  MFEM_ABORT("VTK mesh : inconsistent quadratic mesh!");
-               }
-            }
-         }
-      }
-
-      // Define the 'Nodes' from the 'points' through the 'pts_dof' map
-      for (i = 0; i < np; i++)
-      {
-         dofs.SetSize(1);
-         if ((dofs[0] = pts_dof[i]) != -1)
-         {
-            fes->DofsToVDofs(dofs);
-            for (j = 0; j < dofs.Size(); j++)
-            {
-               (*Nodes)(dofs[j]) = points(3*i+j);
-            }
-         }
-      }
-
-      read_gf = 0;
-   }
-}
+   CreateVTKMesh(points, cell_data, cell_offsets, cell_types, cell_attributes,
+                 curved, read_gf, finalize_topo);
+} // end ReadVTKMesh
 
 void Mesh::ReadNURBSMesh(std::istream &input, int &curved, int &read_gf)
 {
@@ -715,8 +916,8 @@ void Mesh::ReadNURBSMesh(std::istream &input, int &curved, int &read_gf)
       NURBSext->SetCoordsFromPatches(*Nodes);
       own_nodes = 1;
       read_gf = 0;
-      int vd = Nodes->VectorDim();
-      for (int i = 0; i < vd; i++)
+      spaceDim = Nodes->VectorDim();
+      for (int i = 0; i < spaceDim; i++)
       {
          Vector vert_val;
          Nodes->GetNodalValues(vert_val, i+1);

mesh/vtk.cpp

@@ -18,7 +18,137 @@
 namespace mfem
 {
 
-const int vtk_prism_perm[6] = {0, 2, 1, 3, 5, 4};
+const int VTKGeometry::Map[Geometry::NUM_GEOMETRIES] =
+{
+   POINT, SEGMENT, TRIANGLE, SQUARE, TETRAHEDRON, CUBE, PRISM
+};
+
+const int VTKGeometry::QuadraticMap[Geometry::NUM_GEOMETRIES] =
+{
+   POINT, QUADRATIC_SEGMENT, QUADRATIC_TRIANGLE, BIQUADRATIC_SQUARE,
+   QUADRATIC_TETRAHEDRON, TRIQUADRATIC_CUBE, BIQUADRATIC_QUADRATIC_PRISM
+};
+
+const int VTKGeometry::HighOrderMap[Geometry::NUM_GEOMETRIES] =
+{
+   POINT, LAGRANGE_SEGMENT, LAGRANGE_TRIANGLE, LAGRANGE_SQUARE,
+   LAGRANGE_TETRAHEDRON, LAGRANGE_CUBE, LAGRANGE_PRISM
+};
+
+const int VTKGeometry::PrismMap[6] = {0, 2, 1, 3, 5, 4};
+
+const int *VTKGeometry::VertexPermutation[Geometry::NUM_GEOMETRIES] =
+{
+   NULL, NULL, NULL, NULL, NULL, NULL, VTKGeometry::PrismMap
+};
+
+Geometry::Type VTKGeometry::GetMFEMGeometry(int vtk_geom)
+{
+   switch (vtk_geom)
+   {
+      case POINT:
+         return Geometry::POINT;
+      case SEGMENT:
+      case QUADRATIC_SEGMENT:
+      case LAGRANGE_SEGMENT:
+         return Geometry::SEGMENT;
+      case TRIANGLE:
+      case QUADRATIC_TRIANGLE:
+      case LAGRANGE_TRIANGLE:
+         return Geometry::TRIANGLE;
+      case SQUARE:
+      case BIQUADRATIC_SQUARE:
+      case LAGRANGE_SQUARE:
+         return Geometry::SQUARE;
+      case TETRAHEDRON:
+      case QUADRATIC_TETRAHEDRON:
+      case LAGRANGE_TETRAHEDRON:
+         return Geometry::TETRAHEDRON;
+      case CUBE:
+      case TRIQUADRATIC_CUBE:
+      case LAGRANGE_CUBE:
+         return Geometry::CUBE;
+      case PRISM:
+      case BIQUADRATIC_QUADRATIC_PRISM:
+      case LAGRANGE_PRISM:
+         return Geometry::PRISM;
+      default:
+         return Geometry::INVALID;
+   }
+}
+
+bool VTKGeometry::IsLagrange(int vtk_geom)
+{
+   return vtk_geom >= LAGRANGE_SEGMENT && vtk_geom <= LAGRANGE_PRISM;
+}
+
+bool VTKGeometry::IsQuadratic(int vtk_geom)
+{
+   return vtk_geom >= QUADRATIC_SEGMENT
+          && vtk_geom <= BIQUADRATIC_QUADRATIC_PRISM;
+}
+
+int VTKGeometry::GetOrder(int vtk_geom, int npoints)
+{
+   if (IsQuadratic(vtk_geom))
+   {
+      return 2;
+   }
+   else if (IsLagrange(vtk_geom))
+   {
+      switch (vtk_geom)
+      {
+         case LAGRANGE_SEGMENT:
+            return npoints - 1;
+         case LAGRANGE_TRIANGLE:
+            return (std::sqrt(8*npoints + 1) - 3)/2;
+         case LAGRANGE_SQUARE:
+            return std::round(std::sqrt(npoints)) - 1;
+         case LAGRANGE_TETRAHEDRON:
+            switch (npoints)
+            {
+               // Note that for given order, npoints is given by
+               // npoints_order = (order + 1)*(order + 2)*(order + 3)/6,
+               case 4: return 1;
+               case 10: return 2;
+               case 20: return 3;
+               case 35: return 4;
+               case 56: return 5;
+               case 84: return 6;
+               case 120: return 7;
+               case 165: return 8;
+               case 220: return 9;
+               case 286: return 10;
+               default:
+               {
+                  constexpr int max_order = 20;
+                  int order = 11, npoints_order;
+                  for (; order<max_order; ++order)
+                  {
+                     npoints_order = (order + 1)*(order + 2)*(order + 3)/6;
+                     if (npoints_order == npoints) { break; }
+                  }
+                  MFEM_VERIFY(npoints == npoints_order, "");
+                  return order;
+               }
+            }
+         case LAGRANGE_CUBE:
+            return std::round(std::cbrt(npoints)) - 1;
+         case LAGRANGE_PRISM:
+         {
+            const double n = npoints;
+            static const double third = 1.0/3.0;
+            static const double ninth = 1.0/9.0;
+            static const double twentyseventh = 1.0/27.0;
+            const double term =
+               std::cbrt(third*sqrt(third)*sqrt((27.0*n - 2.0)*n) + n
+                         - twentyseventh);
+            return std::round(term + ninth / term - 4*third);
+         }
+      }
+   }
+   return 1;
+}
 
 int BarycentricToVTKTriangle(int *b, int ref)
 {
@@ -379,15 +509,12 @@ void CreateVTKElementConnectivity(Array<int> &con, Geometry::Type geom, int ref)
    {
       int idx = 0;
       int b[4];
-      for (int k=0; k<=ref; k++)
+      for (b[2]=0; b[2]<=ref; b[2]++)
       {
-         for (int j=0; j<=k; j++)
+         for (b[1]=0; b[1]<=ref-b[2]; b[1]++)
          {
-            for (int i=0; i<=j; i++)
+            for (b[0]=0; b[0]<=ref-b[1]-b[2]; b[0]++)
             {
-               b[0] = k-j;
-               b[1] = i;
-               b[2] = j-i;
                b[3] = ref-b[0]-b[1]-b[2];
                con[BarycentricToVTKTetra(b, ref)] = idx++;
             }

mesh/vtk.hpp

@@ -17,7 +17,46 @@
 namespace mfem
 {
 
-// Helpers for writing to the VTK format
+// Helpers for reading and writing VTK format
+
+// VTK element types defined at: https://git.io/JvZLm
+struct VTKGeometry
+{
+   static const int POINT = 1;
+   static const int SEGMENT = 3;
+   static const int TRIANGLE = 5;
+   static const int SQUARE = 9;
+   static const int TETRAHEDRON = 10;
+   static const int CUBE = 12;
+   static const int PRISM = 13;
+
+   static const int QUADRATIC_SEGMENT = 21;
+   static const int QUADRATIC_TRIANGLE = 22;
+   static const int BIQUADRATIC_SQUARE = 28;
+   static const int QUADRATIC_TETRAHEDRON = 24;
+   static const int TRIQUADRATIC_CUBE = 29;
+   static const int QUADRATIC_PRISM = 26;
+   static const int BIQUADRATIC_QUADRATIC_PRISM = 32;
+
+   static const int LAGRANGE_SEGMENT = 68;
+   static const int LAGRANGE_TRIANGLE = 69;
+   static const int LAGRANGE_SQUARE = 70;
+   static const int LAGRANGE_TETRAHEDRON = 71;
+   static const int LAGRANGE_CUBE = 72;
+   static const int LAGRANGE_PRISM = 73;
+
+   static const int PrismMap[6];
+   static const int *VertexPermutation[Geometry::NUM_GEOMETRIES];
+
+   static const int Map[Geometry::NUM_GEOMETRIES];
+   static const int QuadraticMap[Geometry::NUM_GEOMETRIES];
+   static const int HighOrderMap[Geometry::NUM_GEOMETRIES];
+
+   static Geometry::Type GetMFEMGeometry(int vtk_geom);
+   static bool IsLagrange(int vtk_geom);
+   static bool IsQuadratic(int vtk_geom);
+   static int GetOrder(int vtk_geom, int npoints);
+};
 
 enum class VTKFormat
 {
@@ -26,8 +65,6 @@ enum class VTKFormat
    BINARY32
 };
 
-extern const int vtk_prism_perm[6];
-
 /// Create the VTK element connectivity array for a given element geometry and
 /// refinement level. Converts node numbers from MFEM to VTK ordering.
 void CreateVTKElementConnectivity(Array<int> &con, Geometry::Type geom,

tests/unit/CMakeLists.txt

@@ -42,6 +42,7 @@ set(UNIT_TESTS_SRCS
   fem/test_2d_bilininteg.cpp
   fem/test_3d_bilininteg.cpp
   fem/test_assemblediagonalpa.cpp
+  fem/test_blocknonlinearform.cpp
   fem/test_calcshape.cpp
   fem/test_datacollection.cpp
   fem/test_estimator.cpp
@@ -50,6 +51,7 @@ set(UNIT_TESTS_SRCS
   fem/test_intrules.cpp
   fem/test_intruletypes.cpp
   fem/test_inversetransform.cpp
+  fem/test_lexicographic_ordering.cpp
   fem/test_lin_interp.cpp
   fem/test_linear_fes.cpp
   fem/test_operatorjacobismoother.cpp
@@ -59,7 +61,7 @@ set(UNIT_TESTS_SRCS
   fem/test_pa_idinterp.cpp
   fem/test_quadf_coef.cpp
   fem/test_quadraturefunc.cpp
-  fem/test_blocknonlinearform.cpp
+  fem/test_sum_bilin.cpp
   miniapps/test_sedov.cpp
 )
 

tests/unit/fem/test_lexicographic_ordering.cpp

@@ -0,0 +1,64 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+using namespace mfem;
+
+void VerifyOrdering(NodalFiniteElement &el)
+{
+   int order = el.GetOrder();
+   Geometry::Type geom = el.GetGeomType();
+   const Array<int> &p = el.GetLexicographicOrdering();
+
+   GeometryRefiner refiner;
+   refiner.SetType(BasisType::GaussLobatto);
+   RefinedGeometry *ref_geom = refiner.Refine(geom, order);
+
+   double error = 0.0;
+
+   for (int i=0; i<el.GetDof(); ++i)
+   {
+      int pi = (p.Size() > 0) ? p[i] : i;
+      error += std::fabs(el.GetNodes()[pi].x - ref_geom->RefPts[i].x);
+      error += std::fabs(el.GetNodes()[pi].y - ref_geom->RefPts[i].y);
+      error += std::fabs(el.GetNodes()[pi].z - ref_geom->RefPts[i].z);
+   }
+
+   REQUIRE(error == MFEM_Approx(0.0));
+}
+
+template <typename T> void VerifyOrdering(int order)
+{
+   T el(order, BasisType::GaussLobatto);
+   Geometry::Type geom = el.GetGeomType();
+   INFO("order " << order << " " << Geometry::Name[geom]);
+   VerifyOrdering(el);
+}
+
+TEST_CASE("Lexicographic Ordering", "[FiniteElement,Geometry]")
+{
+   auto order = GENERATE(1, 2, 3, 4, 5, 6);
+   VerifyOrdering<H1_SegmentElement>(order);
+   VerifyOrdering<H1_TriangleElement>(order);
+   VerifyOrdering<H1_QuadrilateralElement>(order);
+   VerifyOrdering<H1_TetrahedronElement>(order);
+   VerifyOrdering<H1_HexahedronElement>(order);
+   VerifyOrdering<H1_WedgeElement>(order);
+
+   VerifyOrdering<L2_SegmentElement>(order);
+   VerifyOrdering<L2_TriangleElement>(order);
+   VerifyOrdering<L2_QuadrilateralElement>(order);
+   VerifyOrdering<L2_TetrahedronElement>(order);
+   VerifyOrdering<L2_HexahedronElement>(order);
+   VerifyOrdering<L2_WedgeElement>(order);
+}

tests/unit/fem/test_sum_bilin.cpp

@@ -0,0 +1,231 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "unit_tests.hpp"
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace mfem;
+
+namespace pa_kernels
+{
+
+TEST_CASE("H1 SumIntegrator", "[SumIntegrator][PartialAssembly]")
+{
+   Mesh mesh(1, 1, 1, Element::HEXAHEDRON);
+   H1_FECollection fec(2, mesh.Dimension());
+   FiniteElementSpace fes(&mesh, &fec);
+
+   MassIntegrator integ1;
+   DiffusionIntegrator integ2;
+
+   SumIntegrator integ_sum(true);
+   integ_sum.AddIntegrator(new MassIntegrator);
+   integ_sum.AddIntegrator(new DiffusionIntegrator);
+
+   const FiniteElement &el = *fes.GetFE(0);
+   ElementTransformation &T = *mesh.GetElementTransformation(0);
+   DenseMatrix m1, m_tmp, m2;
+
+   // AssembleElementMatrix
+   integ1.AssembleElementMatrix(el, T, m1);
+   integ2.AssembleElementMatrix(el, T, m_tmp);
+   m1 += m_tmp;
+   integ_sum.AssembleElementMatrix(el, T, m2);
+   m1 -= m2;
+   REQUIRE(m1.MaxMaxNorm() == MFEM_Approx(0.0));
+
+   // AssembleElementMatrix2
+   integ1.AssembleElementMatrix2(el, el, T, m1);
+   integ2.AssembleElementMatrix2(el, el, T, m_tmp);
+   m1 += m_tmp;
+   integ_sum.AssembleElementMatrix2(el, el, T, m2);
+   m1 -= m2;
+   REQUIRE(m1.MaxMaxNorm() == MFEM_Approx(0.0));
+
+   // PA
+   integ1.AssemblePA(fes);
+   integ2.AssemblePA(fes);
+   integ_sum.AssemblePA(fes);
+
+   int n = fes.GetTrueVSize();
+   Vector x(n), y1(n), y2(n);
+   Vector diag1(n), diag_tmp(n), diag2(n);
+   x.Randomize(1);
+
+   // AddMultPA
+   y1 = 0.0;
+   y2 = 0.0;
+   integ1.AddMultPA(x, y1);
+   integ2.AddMultPA(x, y1);
+   integ_sum.AddMultPA(x, y2);
+   y1 -= y2;
+   REQUIRE(y1.Normlinf() == MFEM_Approx(0.0));
+
+   // AddMultTransposePA
+   y1 = 0.0;
+   y2 = 0.0;
+   integ1.AddMultTransposePA(x, y1);
+   integ2.AddMultTransposePA(x, y1);
+   integ_sum.AddMultTransposePA(x, y2);
+   y1 -= y2;
+   REQUIRE(y1.Normlinf() == MFEM_Approx(0.0));
+
+   // AssembleDiagonalPA
+   diag1 = 0.0;
+   diag_tmp = 0.0;
+   diag2 = 0.0;
+   integ1.AssembleDiagonalPA(diag1);
+   integ2.AssembleDiagonalPA(diag_tmp);
+   diag1 += diag_tmp;
+   integ_sum.AssembleDiagonalPA(diag2);
+   diag1 -= diag2;
+   REQUIRE(diag1.Normlinf() == MFEM_Approx(0.0));
+
+   // MF
+#ifdef MFEM_USE_CEED
+   if (DeviceCanUseCeed())
+   {
+      integ1.AssembleMF(fes);
+      integ2.AssembleMF(fes);
+      integ_sum.AssembleMF(fes);
+
+      // AddMultMF
+      y1 = 0.0;
+      y2 = 0.0;
+      integ1.AddMultMF(x, y1);
+      integ2.AddMultMF(x, y1);
+      integ_sum.AddMultMF(x, y2);
+      y1 -= y2;
+      REQUIRE(y1.Normlinf() == MFEM_Approx(0.0));
+
+      // AddMultTransposeMF
+      y1 = 0.0;
+      y2 = 0.0;
+      integ1.AddMultTransposeMF(x, y1);
+      integ2.AddMultTransposeMF(x, y1);
+      integ_sum.AddMultTransposeMF(x, y2);
+      y1 -= y2;
+      REQUIRE(y1.Normlinf() == MFEM_Approx(0.0));
+
+      // AssembleDiagonalMF
+      integ1.AssembleDiagonalMF(diag1);
+      integ2.AssembleDiagonalMF(diag_tmp);
+      diag1 += diag_tmp;
+      integ_sum.AssembleDiagonalMF(diag2);
+      diag1 -= diag2;
+      REQUIRE(diag1.Normlinf() == MFEM_Approx(0.0));
+   }
+#endif
+}
+
+TEST_CASE("DG SumIntegrator", "[SumIntegrator][PartialAssembly]")
+{
+   Mesh mesh(2, 1, 1, Element::HEXAHEDRON);
+   DG_FECollection fec(2, mesh.Dimension(), BasisType::GaussLobatto);
+   FiniteElementSpace fes(&mesh, &fec);
+
+   Vector v(mesh.Dimension());
+   v = 1.0;
+   VectorConstantCoefficient v_coeff(v);
+
+   DGTraceIntegrator integ1(v_coeff, 1.0, 2.0);
+   DGTraceIntegrator integ2(v_coeff, 3.0, 4.0);
+
+   SumIntegrator integ_sum(true);
+   integ_sum.AddIntegrator(new DGTraceIntegrator(v_coeff, 1.0, 2.0));
+   integ_sum.AddIntegrator(new DGTraceIntegrator(v_coeff, 3.0, 4.0));
+
+   DenseMatrix m1, m_tmp, m2;
+
+   // AssembleFaceMatrix
+   int nfaces = mesh.GetNumFaces();
+   for (int i = 0; i < nfaces; i++)
+   {
+      FaceElementTransformations *tr = mesh.GetFaceElementTransformations(i);
+      const FiniteElement &el0 = *fes.GetFE(tr->Elem1No);
+      const FiniteElement &el1 = (tr->Elem2No >= 0) ? *fes.GetFE(tr->Elem2No) : el0;
+      integ1.AssembleFaceMatrix(el0, el1, *tr, m1);
+      integ2.AssembleFaceMatrix(el0, el1, *tr, m_tmp);
+      m1 += m_tmp;
+      integ_sum.AssembleFaceMatrix(el0, el1, *tr, m2);
+      m1 -= m2;
+      REQUIRE(m1.MaxMaxNorm() == MFEM_Approx(0.0));
+   }
+
+   // PA interior
+   integ1.AssemblePAInteriorFaces(fes);
+   integ2.AssemblePAInteriorFaces(fes);
+   integ_sum.AssemblePAInteriorFaces(fes);
+
+   const Operator *R_int = fes.GetFaceRestriction(
+                              ElementDofOrdering::LEXICOGRAPHIC,
+                              FaceType::Interior);
+
+   int n_int = R_int->Height();
+   Vector x(n_int), y1(n_int), y2(n_int);
+   x.Randomize(1);
+
+   // AddMultPA
+   y1 = 0.0;
+   y2 = 0.0;
+   integ1.AddMultPA(x, y1);
+   integ2.AddMultPA(x, y1);
+   integ_sum.AddMultPA(x, y2);
+   y1 -= y2;
+   REQUIRE(y1.Normlinf() == MFEM_Approx(0.0));
+
+   // AddMultTransposePA
+   y1 = 0.0;
+   y2 = 0.0;
+   integ1.AddMultTransposePA(x, y1);
+   integ2.AddMultTransposePA(x, y1);
+   integ_sum.AddMultTransposePA(x, y2);
+   y1 -= y2;
+   REQUIRE(y1.Normlinf() == MFEM_Approx(0.0));
+
+   // PA boundary
+   integ1.AssemblePABoundaryFaces(fes);
+   integ2.AssemblePABoundaryFaces(fes);
+   integ_sum.AssemblePABoundaryFaces(fes);
+
+   const Operator *R_bdr = fes.GetFaceRestriction(
+                              ElementDofOrdering::LEXICOGRAPHIC,
+                              FaceType::Boundary,
+                              L2FaceValues::DoubleValued);
+
+   int n_bdr = R_bdr->Height();
+   x.SetSize(n_bdr);
+   y1.SetSize(n_bdr);
+   y2.SetSize(n_bdr);
+   x.Randomize(1);
+
+   // AddMultPA
+   y1 = 0.0;
+   y2 = 0.0;
+   integ1.AddMultPA(x, y1);
+   integ2.AddMultPA(x, y1);
+   integ_sum.AddMultPA(x, y2);
+   y1 -= y2;
+   REQUIRE(y1.Normlinf() == MFEM_Approx(0.0));
+
+   // AddMultTransposePA
+   y1 = 0.0;
+   y2 = 0.0;
+   integ1.AddMultTransposePA(x, y1);
+   integ2.AddMultTransposePA(x, y1);
+   integ_sum.AddMultTransposePA(x, y2);
+   y1 -= y2;
+   REQUIRE(y1.Normlinf() == MFEM_Approx(0.0));
+}
+
+} // namespace pa_kernels