Updated explanation from PR

CHANGELOG

@@ -11,9 +11,6 @@
 Version 4.2.1 (development)
 ===========================
 
-- Added matrix-free GPU-enabled implementations of GradientInterpolator and
-  IdentityInterpolator.
-
 - Added interface to MUMPS direct solver. Its usage is demonstrated in ex25p.
   See http://mumps.enseeiht.fr/ for more details. Supported versions >= 5.1.1.
 
@@ -35,12 +32,6 @@ 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.
 
@@ -58,11 +49,6 @@ 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,5 +29,4 @@ 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,9 +819,7 @@ RECURSIVE              = NO
 # run.
 
 EXCLUDE                = @MFEM_SOURCE_DIR@/config/_config.hpp \
-                         @MFEM_SOURCE_DIR@/config/get_hypre_version.cpp \
-                         @MFEM_SOURCE_DIR@/general/tinyxml2.h \
-                         @MFEM_SOURCE_DIR@/general/tinyxml2.cpp
+                         @MFEM_SOURCE_DIR@/config/get_hypre_version.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,7 +34,6 @@ list(APPEND ALL_EXE_SRCS
   ex25.cpp
   ex26.cpp
   ex27.cpp
-  ex91.cpp
   )
 
 if (MFEM_USE_MPI)

examples/ex16p.cpp

@@ -65,7 +65,7 @@ protected:
    HypreSmoother M_prec; // Preconditioner for the mass matrix M
 
    CGSolver T_solver;    // Implicit solver for T = M + dt K
-   HypreSmoother T_prec; // Preconditioner for the implicit solver
+   HypreBoomerAMG T_prec; // Preconditioner for the implicit solver
 
    double alpha, kappa;
 
@@ -353,7 +353,6 @@ int main(int argc, char *argv[])
          }
 #endif
       }
-      oper.SetParameters(u);
    }
 
 #ifdef MFEM_USE_ADIOS2
@@ -414,6 +413,7 @@ ConductionOperator::ConductionOperator(ParFiniteElementSpace &f, double al,
    T_solver.SetPrintLevel(0);
    T_solver.SetPreconditioner(T_prec);
 
+   T_prec.SetPrintLevel(0);
    SetParameters(u);
 }
 
@@ -430,19 +430,53 @@ void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
 void ConductionOperator::ImplicitSolve(const double dt,
                                        const Vector &u, Vector &du_dt)
 {
-   // Solve the equation:
-   //    du_dt = M^{-1}*[-K(u + dt*du_dt)]
-   // for du_dt
-   if (!T)
-   {
-      T = Add(1.0, Mmat, dt, Kmat);
-      current_dt = dt;
+   // Here we use Picard iterations to solve a nonlinear equation
+   // for the Runge-Kutta stage vector k,
+   //
+   //    M*k = N(u+dt*k)                     (1)
+   //
+   // for nonlinear operator N. We assume N can be written as
+   //
+   //    N(u+dt*k) := L[u+dt*k](u+dt*k) + f(t)
+   //
+   // where L is a matrix-valued operator evaluated at u+dt*k and f(t)
+   // a (potentially zero) time-dependent forcing vector. (1) can be
+   // rewritten as a fixed-point equation
+   //
+   //    x = (M - dt*L[x])^{-1} (Mu + f)     (2)
+   //
+   // where x := u + dt*k, which can be solved using a Picard iteration,
+   // where a function G(x) = x is solved via iteraitons x_{k+1} = G(x_k).
+
+   double tol = 1e-6;
+   int maxiter = 100;
+
+   // Right-hand side for nonlinear iteration
+   Mmat.Mult(u, z);  // Add forcing function if one exists
+   du_dt = u;        // Set u as initial guess for x (2)
+   Vector temp(u);   // Vector to measure error
+   temp = u;
+   double error = 1;
+   int iter = 0;
+   while (error > tol) {
+      iter ++;
+      this->SetParameters(du_dt);   // Update nonlinear operator L[x]
+      T = Add(1.0, Mmat, dt, Kmat); // Form matrix (M - dt*L[x])
       T_solver.SetOperator(*T);
+      T_solver.Mult(z, du_dt);      // Apply (M - dt*L[x])^{-1}
+      temp -= du_dt;                // Measure error
+      error = std::sqrt(InnerProduct(MPI_COMM_WORLD, temp, temp));
+      temp = du_dt;
+      if (iter >= maxiter) {
+         mfem_warning("Nonlinear iteration did not converge!");
+         break;
+      }
    }
-   MFEM_VERIFY(dt == current_dt, ""); // SDIRK methods use the same dt
-   Kmat.Mult(u, z);
-   z.Neg();
-   T_solver.Mult(z, du_dt);
+
+   // Above we solved for x = u + dt*k, where k is the desired update
+   // Map du_dt -> k.
+   du_dt -= u;
+   du_dt /= dt;
 }
 
 void ConductionOperator::SetParameters(const Vector &u)
@@ -483,4 +517,4 @@ double InitialTemperature(const Vector &x)
    {
       return 1.0;
    }
-}
+}

examples/ex1p.cpp

@@ -231,10 +231,17 @@ int main(int argc, char *argv[])
    else
    {
       prec = new HypreBoomerAMG;
+      prec->SetOperator(*A);
    }
-   CGSolver cg(MPI_COMM_WORLD);
+   // CGSolver cg(MPI_COMM_WORLD);
+   AndersonAcceleration cg(MPI_COMM_WORLD);
+   cg.SetKDim(10);
+   cg.SetRestart(true);      // WORKS
+   cg.SetAAStart(0);          // WORKS
+   cg.SetWeight(1.0);         // Not robust but seems to work
+
    cg.SetRelTol(1e-12);
-   cg.SetMaxIter(2000);
+   cg.SetMaxIter(50);
    cg.SetPrintLevel(1);
    if (prec) { cg.SetPreconditioner(*prec); }
    cg.SetOperator(*A);

examples/ex9.cpp

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

examples/ex91.cpp

@@ -1,672 +0,0 @@
-#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,66 +64,6 @@ 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:
@@ -131,37 +71,24 @@ private:
    SparseMatrix M_diag;
    HypreParMatrix *A;
    GMRESSolver linear_solver;
-   Solver *prec;
+   BlockILU prec;
    double dt;
 public:
-   DG_Solver(HypreParMatrix &M_, HypreParMatrix &K_, const FiniteElementSpace &fes,
-             PrecType prec_type)
+   DG_Solver(HypreParMatrix &M_, HypreParMatrix &K_, const FiniteElementSpace &fes)
       : 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);
    }
@@ -194,12 +121,10 @@ 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
@@ -217,8 +142,7 @@ private:
    mutable Vector z;
 
 public:
-   FE_Evolution(ParBilinearForm &_M, ParBilinearForm &_K, const Vector &_b,
-                PrecType prec_type);
+   FE_Evolution(ParBilinearForm &_M, ParBilinearForm &_K, const Vector &_b);
 
    virtual void Mult(const Vector &x, Vector &y) const;
    virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);
@@ -230,9 +154,10 @@ public:
 int main(int argc, char *argv[])
 {
    // 1. Initialize MPI.
-   MPI_Session mpi;
-   int num_procs = mpi.WorldSize();
-   int myid = mpi.WorldRank();
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
 
    // 2. Parse command-line options.
    problem = 0;
@@ -253,11 +178,7 @@ 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);
 
@@ -291,8 +212,6 @@ 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.");
@@ -313,19 +232,20 @@ int main(int argc, char *argv[])
    args.Parse();
    if (!args.Good())
    {
-      if (mpi.Root())
+      if (myid == 0)
       {
          args.PrintUsage(cout);
       }
+      MPI_Finalize();
       return 1;
    }
-   if (mpi.Root())
+   if (myid == 0)
    {
       args.PrintOptions(cout);
    }
 
    Device device(device_config);
-   if (mpi.Root()) { device.Print(); }
+   if (myid == 0) { device.Print(); }
 
    // 3. Read the serial mesh from the given mesh file on all processors. We can
    //    handle geometrically periodic meshes in this code.
@@ -352,11 +272,12 @@ int main(int argc, char *argv[])
       case 23: ode_solver = new SDIRK23Solver; break;
       case 24: ode_solver = new SDIRK34Solver; break;
       default:
-         if (mpi.Root())
+         if (myid == 0)
          {
             cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
          }
          delete mesh;
+         MPI_Finalize();
          return 3;
    }
 
@@ -390,7 +311,7 @@ int main(int argc, char *argv[])
    ParFiniteElementSpace *fes = new ParFiniteElementSpace(pmesh, &fec);
 
    HYPRE_Int global_vSize = fes->GlobalTrueVSize();
-   if (mpi.Root())
+   if (myid == 0)
    {
       cout << "Number of unknowns: " << global_vSize << endl;
    }
@@ -421,16 +342,15 @@ int main(int argc, char *argv[])
    }
 
    m->AddDomainIntegrator(new MassIntegrator);
-   constexpr double alpha = -1.0;
-   k->AddDomainIntegrator(new ConvectionIntegrator(velocity, alpha));
+   k->AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
    k->AddInteriorFaceIntegrator(
-      new NonconservativeDGTraceIntegrator(velocity, alpha));
+      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
    k->AddBdrFaceIntegrator(
-      new NonconservativeDGTraceIntegrator(velocity, alpha));
+      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
 
    ParLinearForm *b = new ParLinearForm(fes);
    b->AddBdrFaceIntegrator(
-      new BoundaryFlowIntegrator(inflow, velocity, alpha));
+      new BoundaryFlowIntegrator(inflow, velocity, -1.0, -0.5));
 
    int skip_zeros = 0;
    m->Assemble();
@@ -531,11 +451,11 @@ int main(int argc, char *argv[])
       sout.open(vishost, visport);
       if (!sout)
       {
-         if (mpi.Root())
+         if (myid == 0)
             cout << "Unable to connect to GLVis server at "
                  << vishost << ':' << visport << endl;
          visualization = false;
-         if (mpi.Root())
+         if (myid == 0)
          {
             cout << "GLVis visualization disabled.\n";
          }
@@ -547,7 +467,7 @@ int main(int argc, char *argv[])
          sout << "solution\n" << *pmesh << *u;
          sout << "pause\n";
          sout << flush;
-         if (mpi.Root())
+         if (myid == 0)
             cout << "GLVis visualization paused."
                  << " Press space (in the GLVis window) to resume it.\n";
       }
@@ -556,7 +476,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, prec_type);
+   FE_Evolution adv(*m, *k, *B);
 
    double t = 0.0;
    adv.SetTime(t);
@@ -573,7 +493,7 @@ int main(int argc, char *argv[])
 
       if (done || ti % vis_steps == 0)
       {
-         if (mpi.Root())
+         if (myid == 0)
          {
             cout << "time step: " << ti << ", time: " << t << endl;
          }
@@ -644,14 +564,16 @@ 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, PrecType prec_type)
-   : TimeDependentOperator(_M.Height()), b(_b),
+                           const Vector &_b)
+   : TimeDependentOperator(_M.Height()),
+     b(_b),
      M_solver(_M.ParFESpace()->GetComm()),
      z(_M.Height())
 {
@@ -676,7 +598,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(), prec_type);
+      dg_solver = new DG_Solver(M_mat, K_mat, *_M.FESpace());
    }
    else
    {
@@ -692,10 +614,6 @@ 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,7 +66,6 @@ set(SRCS
   tmop_tools.cpp
   gslib.cpp
   transfer.cpp
-  prmnonlinearform.cpp
   )
 
 set(HDRS
@@ -113,7 +112,6 @@ set(HDRS
   tmop_tools.hpp
   gslib.hpp
   transfer.hpp
-  prmnonlinearform.hpp
   )
 
 if (MFEM_USE_SIDRE)
@@ -138,8 +136,7 @@ if (MFEM_USE_MPI)
     pgridfunc.cpp
     plinearform.cpp
     pnonlinearform.cpp
-    prestriction.cpp
-    pprmnonlinearform.cpp)
+    prestriction.cpp)
   # If this list (HDRS -> HEADERS) is used for install, we probably want the
   # headers added all the time.
   list(APPEND HDRS
@@ -148,8 +145,7 @@ if (MFEM_USE_MPI)
     pgridfunc.hpp
     plinearform.hpp
     pnonlinearform.hpp
-    prestriction.hpp
-    pprmnonlinearform.hpp)
+    prestriction.hpp)
 endif()
 
 convert_filenames_to_full_paths(SRCS)

fem/bilinearform.cpp

@@ -1770,41 +1770,9 @@ MixedBilinearForm::~MixedBilinearForm()
    delete ext;
 }
 
-void DiscreteLinearOperator::SetAssemblyLevel(AssemblyLevel assembly_level)
-{
-   if (ext)
-   {
-      MFEM_ABORT("the assembly level has already been set!");
-   }
-   assembly = assembly_level;
-   switch (assembly)
-   {
-      case AssemblyLevel::LEGACYFULL:
-      case AssemblyLevel::FULL:
-         // Use the original implementation for now
-         break;
-      case AssemblyLevel::ELEMENT:
-         mfem_error("Element assembly not supported yet... stay tuned!");
-         break;
-      case AssemblyLevel::PARTIAL:
-         ext = new PADiscreteLinearOperatorExtension(this);
-         break;
-      case AssemblyLevel::NONE:
-         mfem_error("Matrix-free action not supported yet... stay tuned!");
-         break;
-      default:
-         mfem_error("Unknown assembly level");
-   }
-}
 
 void DiscreteLinearOperator::Assemble(int skip_zeros)
 {
-   if (ext)
-   {
-      ext->Assemble();
-      return;
-   }
-
    Array<int> dom_vdofs, ran_vdofs;
    ElementTransformation *T;
    const FiniteElement *dom_fe, *ran_fe;

fem/bilinearform.hpp

@@ -376,13 +376,6 @@ public:
    /// Get the output finite element space prolongation matrix
    virtual const Operator *GetOutputProlongation() const
    { return GetProlongation(); }
-   /** @brief Returns the output fe space restriction matrix, transposed
-
-       Logically, this is the transpose of GetOutputRestriction, but in
-       practice it is convenient to have it in transposed form for
-       construction of RAP operators in matrix-free methods. */
-   virtual const Operator *GetOutputRestrictionTranspose() const
-   { return GetOutputProlongation(); }
    /// Get the output finite element space restriction matrix
    virtual const Operator *GetOutputRestriction() const
    { return GetRestriction(); }
@@ -854,9 +847,9 @@ public:
 
       This returns the same operator as FormRectangularLinearSystem(), but does
       without the transformations of the right-hand side. */
-   virtual void FormRectangularSystemMatrix(const Array<int> &trial_tdof_list,
-                                            const Array<int> &test_tdof_list,
-                                            OperatorHandle &A);
+   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,11 +876,10 @@ 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. */
-   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);
+   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(.).
@@ -985,18 +977,9 @@ public:
    /// Access all interpolators added with AddDomainInterpolator().
    Array<BilinearFormIntegrator*> *GetDI() { return &dbfi; }
 
-   /// Set the desired assembly level. The default is AssemblyLevel::FULL.
-   /** This method must be called before assembly. */
-   void SetAssemblyLevel(AssemblyLevel assembly_level);
-
    /** @brief Construct the internal matrix representation of the discrete
        linear operator. */
    virtual void Assemble(int skip_zeros = 1);
-
-   /** @brief Get the output finite element space restriction matrix in
-       transposed form. */
-   virtual const Operator *GetOutputRestrictionTranspose() const
-   { return test_fes->GetRestrictionTransposeOperator(); }
 };
 
 }

fem/bilinearform_ext.cpp

@@ -1021,6 +1021,7 @@ void PAMixedBilinearFormExtension::Update()
       localTrial.UseDevice(true);
       localTrial.SetSize(elem_restrict_trial->Height(),
                          Device::GetMemoryType());
+
    }
    if (elem_restrict_test)
    {
@@ -1220,134 +1221,4 @@ void PAMixedBilinearFormExtension::AssembleDiagonal_ADAt(const Vector &D,
    }
 }
 
-PADiscreteLinearOperatorExtension::PADiscreteLinearOperatorExtension(
-   DiscreteLinearOperator *linop) :
-   PAMixedBilinearFormExtension(linop)
-{
-}
-
-const
-Operator *PADiscreteLinearOperatorExtension::GetOutputRestrictionTranspose()
-const
-{
-   return a->GetOutputRestrictionTranspose();
-}
-
-void PADiscreteLinearOperatorExtension::Assemble()
-{
-   Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
-   const int integratorCount = integrators.Size();
-   for (int i = 0; i < integratorCount; ++i)
-   {
-      integrators[i]->AssemblePA(*trialFes, *testFes);
-   }
-
-   test_multiplicity.UseDevice(true);
-   test_multiplicity.SetSize(elem_restrict_test->Width()); // l-vector
-   Vector ones(elem_restrict_test->Height()); // e-vector
-   ones = 1.0;
-
-   const ElementRestriction* elem_restrict =
-      dynamic_cast<const ElementRestriction*>(elem_restrict_test);
-   if (elem_restrict)
-   {
-      elem_restrict->MultTransposeUnsigned(ones, test_multiplicity);
-   }
-   else
-   {
-      mfem_error("A real ElementRestriction is required in this setting!");
-   }
-
-   auto tm = test_multiplicity.ReadWrite();
-   MFEM_FORALL(i, test_multiplicity.Size(),
-   {
-      tm[i] = 1.0 / tm[i];
-   });
-}
-
-void PADiscreteLinearOperatorExtension::AddMult(
-   const Vector &x, Vector &y, const double c) const
-{
-   Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
-   const int iSz = integrators.Size();
-
-   // * G operation
-   SetupMultInputs(elem_restrict_trial, x, localTrial,
-                   elem_restrict_test, y, localTest, c);
-
-   // * B^TDB operation
-   for (int i = 0; i < iSz; ++i)
-   {
-      integrators[i]->AddMultPA(localTrial, localTest);
-   }
-
-   // do a kind of "set" rather than "add" in the below
-   // operation as compared to the BilinearForm case
-   // * G^T operation (kind of...)
-   const ElementRestriction* elem_restrict =
-      dynamic_cast<const ElementRestriction*>(elem_restrict_test);
-   if (elem_restrict)
-   {
-      tempY.SetSize(y.Size());
-      elem_restrict->MultLeftInverse(localTest, tempY);
-      y += tempY;
-   }
-   else
-   {
-      mfem_error("In this setting you need a real ElementRestriction!");
-   }
-}
-
-void PADiscreteLinearOperatorExtension::AddMultTranspose(
-   const Vector &x, Vector &y, const double c) const
-{
-   Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
-   const int iSz = integrators.Size();
-
-   // do a kind of "set" rather than "add" in the below
-   // operation as compared to the BilinearForm case
-   // * G operation (kinda)
-   Vector xscaled(x);
-   MFEM_VERIFY(x.Size() == test_multiplicity.Size(), "Input vector of wrong size");
-   auto xs = xscaled.ReadWrite();
-   auto tm = test_multiplicity.Read();
-   MFEM_FORALL(i, x.Size(),
-   {
-      xs[i] *= tm[i];
-   });
-   SetupMultInputs(elem_restrict_test, xscaled, localTest,
-                   elem_restrict_trial, y, localTrial, c);
-
-   // * B^TD^TB operation
-   for (int i = 0; i < iSz; ++i)
-   {
-      integrators[i]->AddMultTransposePA(localTest, localTrial);
-   }
-
-   // * G^T operation
-   if (elem_restrict_trial)
-   {
-      tempY.SetSize(y.Size());
-      elem_restrict_trial->MultTranspose(localTrial, tempY);
-      y += tempY;
-   }
-   else
-   {
-      mfem_error("Trial ElementRestriction not defined");
-   }
-}
-
-void PADiscreteLinearOperatorExtension::FormRectangularSystemOperator(
-   const Array<int>& ess1, const Array<int>& ess2, OperatorHandle &A)
-{
-   const Operator *Pi = this->GetProlongation();
-   const Operator *RoT = this->GetOutputRestrictionTranspose();
-   Operator *rap = SetupRAP(Pi, RoT);
-
-   RectangularConstrainedOperator *Arco
-      = new RectangularConstrainedOperator(rap, ess1, ess2, rap != this);
-
-   A.Reset(Arco);
-}
-
 } // namespace mfem

fem/bilinearform_ext.hpp

@@ -21,7 +21,6 @@ namespace mfem
 
 class BilinearForm;
 class MixedBilinearForm;
-class DiscreteLinearOperator;
 
 /// Class extending the BilinearForm class to support different AssemblyLevels.
 /**  FA - Full Assembly
@@ -213,7 +212,7 @@ protected:
    mutable Vector localTrial, localTest, tempY;
    const Operator *elem_restrict_trial; // Not owned
    const Operator *elem_restrict_test;  // Not owned
-
+private:
    /// Helper function to set up inputs/outputs for Mult or MultTranspose
    void SetupMultInputs(const Operator *elem_restrict_x,
                         const Vector &x, Vector &localX,
@@ -259,35 +258,6 @@ public:
    void Update();
 };
 
-
-/**
-   @brief Partial assembly extension for DiscreteLinearOperator
-
-   This acts very much like PAMixedBilinearFormExtension, but its
-   FormRectangularSystemOperator implementation emulates 'Set' rather than
-   'Add' in the assembly case.
-*/
-class PADiscreteLinearOperatorExtension : public PAMixedBilinearFormExtension
-{
-public:
-   PADiscreteLinearOperatorExtension(DiscreteLinearOperator *linop);
-
-   /// Partial assembly of all internal integrators
-   void Assemble();
-
-   void AddMult(const Vector &x, Vector &y, const double c) const;
-
-   void AddMultTranspose(const Vector &x, Vector &y, const double c=1.0) const;
-
-   void FormRectangularSystemOperator(const Array<int>&, const Array<int>&,
-                                      OperatorHandle& A);
-
-   const Operator * GetOutputRestrictionTranspose() const;
-
-private:
-   Vector test_multiplicity;
-};
-
 }
 
 #endif

fem/bilininteg.cpp

@@ -22,14 +22,14 @@ namespace mfem
 
 void BilinearFormIntegrator::AssemblePA(const FiniteElementSpace&)
 {
-   mfem_error ("BilinearFormIntegrator::AssemblePA(fes)\n"
+   mfem_error ("BilinearFormIntegrator::AssemblePA(...)\n"
                "   is not implemented for this class.");
 }
 
 void BilinearFormIntegrator::AssemblePA(const FiniteElementSpace&,
                                         const FiniteElementSpace&)
 {
-   mfem_error ("BilinearFormIntegrator::AssemblePA(fes, fes)\n"
+   mfem_error ("BilinearFormIntegrator::AssemblePA(...)\n"
                "   is not implemented for this class.");
 }
 
@@ -92,7 +92,7 @@ void BilinearFormIntegrator::AddMultPA(const Vector &, Vector &) const
 
 void BilinearFormIntegrator::AddMultTransposePA(const Vector &, Vector &) const
 {
-   mfem_error ("BilinearFormIntegrator::AddMultTransposePA(...)\n"
+   mfem_error ("BilinearFormIntegrator::MultAssembledTranspose(...)\n"
                "   is not implemented for this class.");
 }
 
@@ -229,159 +229,6 @@ 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;
-   mutable DenseMatrix elem_mat;
+   DenseMatrix elem_mat;
    Array<BilinearFormIntegrator*> integrators;
 
 public:
@@ -367,55 +367,6 @@ 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();
 };
@@ -1893,10 +1844,8 @@ protected:
 };
 
 /** Class for integrating the bilinear form a(u,v) := (Q grad u, v) where Q is a
-    scalar coefficient, and v is a vector with components v_i in the same (H1) space
-    as u.
-
-    See also MixedVectorGradientIntegrator when v is in H(curl). */
+    scalar coefficient, and v is a vector with components v_i in the same space
+    as u. */
 class GradientIntegrator : public BilinearFormIntegrator
 {
 protected:
@@ -2045,8 +1994,6 @@ 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);
 };
@@ -2110,8 +2057,6 @@ 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);
@@ -2173,17 +2118,6 @@ 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
 {
@@ -2748,15 +2682,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, see also
-    NonconservativeDGTraceIntegrator). 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). 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 ConservativeConvectionIntegrator (with coefficient u, and
-    parameter alpha = -1) plus the DGTraceIntegrator (with coefficient u, and
+    consists of the transpose of the ConvectionIntegrator (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
@@ -2775,17 +2709,13 @@ private:
    Vector shape1, shape2;
 
 public:
-   /// Construct integrator with rho = 1, b = 0.5*a.
-   DGTraceIntegrator(VectorCoefficient &u_, double a)
-   { rho = NULL; u = &u_; alpha = a; beta = 0.5*a; }
-
    /// Construct integrator with rho = 1.
-   DGTraceIntegrator(VectorCoefficient &u_, double a, double b)
-   { rho = NULL; u = &u_; alpha = a; beta = b; }
+   DGTraceIntegrator(VectorCoefficient &_u, double a, double b)
+   { rho = NULL; u = &_u; alpha = a; beta = b; }
 
-   DGTraceIntegrator(Coefficient &_rho, VectorCoefficient &u_,
+   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,
@@ -2819,30 +2749,6 @@ 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} >
@@ -3083,36 +2989,11 @@ class DiscreteInterpolator : public BilinearFormIntegrator { };
 class GradientInterpolator : public DiscreteInterpolator
 {
 public:
-   GradientInterpolator() : dofquad_fe(NULL) { }
-   virtual ~GradientInterpolator() { delete dofquad_fe; }
-
    virtual void AssembleElementMatrix2(const FiniteElement &h1_fe,
                                        const FiniteElement &nd_fe,
                                        ElementTransformation &Trans,
                                        DenseMatrix &elmat)
    { nd_fe.ProjectGrad(h1_fe, Trans, elmat); }
-
-   using BilinearFormIntegrator::AssemblePA;
-
-   /** @brief Setup method for PA data.
-
-       @param[in] trial_fes   H1 Lagrange space
-       @param[in] test_fes    H(curl) Nedelec space
-    */
-   virtual void AssemblePA(const FiniteElementSpace &trial_fes,
-                           const FiniteElementSpace &test_fes);
-
-   virtual void AddMultPA(const Vector &x, Vector &y) const;
-   virtual void AddMultTransposePA(const Vector &x, Vector &y) const;
-
-private:
-   /// 1D finite element that generates and owns the 1D DofToQuad maps below
-   FiniteElement * dofquad_fe;
-
-   bool B_id; // is the B basis operator (maps_C_C) the identity?
-   const DofToQuad *maps_C_C; // one-d map with Lobatto rows, Lobatto columns
-   const DofToQuad *maps_O_C; // one-d map with Legendre rows, Lobatto columns
-   int dim, ne, o_dofs1D, c_dofs1D;
 };
 
 
@@ -3127,24 +3008,6 @@ public:
                                        ElementTransformation &Trans,
                                        DenseMatrix &elmat)
    { ran_fe.Project(dom_fe, Trans, elmat); }
-
-   using BilinearFormIntegrator::AssemblePA;
-
-   virtual void AssemblePA(const FiniteElementSpace &trial_fes,
-                           const FiniteElementSpace &test_fes);
-
-   virtual void AddMultPA(const Vector &x, Vector &y) const;
-   virtual void AddMultTransposePA(const Vector &x, Vector &y) const;
-
-private:
-   /// 1D finite element that generates and owns the 1D DofToQuad maps below
-   FiniteElement * dofquad_fe;
-
-   const DofToQuad *maps_C_C; // one-d map with Lobatto rows, Lobatto columns
-   const DofToQuad *maps_O_C; // one-d map with Legendre rows, Lobatto columns
-   int dim, ne, o_dofs1D, c_dofs1D;
-
-   Vector pa_data;
 };
 
 

fem/bilininteg_diffusion_pa.cpp

@@ -1904,17 +1904,4 @@ 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_gradient.cpp

@@ -334,7 +334,7 @@ static void PAGradientApplyTranspose2D(const int NE,
                                        const int q1d = 0)
 {
    // TODO
-   MFEM_ASSERT(false, "PAGradientApplyTranspose2D not implemented.");
+   MFEM_ASSERT(false, "GradientPAApplyTranspose 3D not implemented.");
 }
 
 // PA Gradient Apply 3D kernel

fem/bilininteg_mass_pa.cpp

@@ -1226,10 +1226,4 @@ 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,10 +7797,7 @@ 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)
-{
-   lex_ordering = dof_map;
-}
+     TensorBasisElement(dims, p, VerifyNodal(btype), dmtype) { }
 
 
 PositiveTensorFiniteElement::PositiveTensorFiniteElement(
@@ -8484,33 +8481,23 @@ 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]);
    }
 
@@ -8519,7 +8506,6 @@ 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);
       }
 
@@ -8653,56 +8639,36 @@ 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]);
    }
 
@@ -8710,28 +8676,24 @@ 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]);
       }
@@ -8741,7 +8703,6 @@ 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);
          }
@@ -9286,22 +9247,7 @@ 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;
@@ -9310,19 +9256,9 @@ 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;
@@ -9332,26 +9268,21 @@ 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
-   k=0;
+   int 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;
@@ -9360,10 +9291,6 @@ 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;
@@ -9375,7 +9302,6 @@ H1_WedgeElement::H1_WedgeElement(const int p,
          k++;
       }
    }
-   o += 3*nq;
 
    // Interior DoFs
    int m=0;
@@ -9384,9 +9310,8 @@ H1_WedgeElement::H1_WedgeElement(const int p,
       int l=0;
       for (int j=1; j<p; j++)
       {
-         for (int i=1; i+j<p; i++)
+         for (int i=1; i<j; 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,7 +698,6 @@ public:
 class NodalFiniteElement : public ScalarFiniteElement
 {
 protected:
-   Array<int> lex_ordering;
    void ProjectCurl_2D(const FiniteElement &fe,
                        ElementTransformation &Trans,
                        DenseMatrix &curl) const;
@@ -747,29 +746,6 @@ 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,16 +1854,19 @@ const int *H1_FECollection::GetDofMap(Geometry::Type GeomType) const
 {
    const int *dof_map = NULL;
    const FiniteElement *fe = H1_Elements[GeomType];
-   const NodalFiniteElement *nodal_fe =
-      dynamic_cast<const NodalFiniteElement*>(fe);
-   if (nodal_fe)
+   switch (GeomType)
    {
-      dof_map = nodal_fe->GetLexicographicOrdering().GetData();
-   }
-   else
-   {
-      MFEM_ABORT("Geometry type " << Geometry::Name[GeomType] << " is not "
-                 "implemented");
+      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.
    }
    return dof_map;
 }

fem/fem.hpp

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

fem/fespace.cpp

@@ -2608,9 +2608,9 @@ const Operator &GridTransfer::MakeTrueOperator(
    else // Parallel() == true
    {
 #ifdef MFEM_USE_MPI
+      const SparseMatrix *out_R = fes_out.GetRestrictionMatrix();
       if (oper_type == Operator::Hypre_ParCSR)
       {
-         const SparseMatrix *out_R = fes_out.GetRestrictionMatrix();
          const ParFiniteElementSpace *pfes_in =
             dynamic_cast<const ParFiniteElementSpace *>(&fes_in);
          const ParFiniteElementSpace *pfes_out =
@@ -2638,7 +2638,6 @@ const Operator &GridTransfer::MakeTrueOperator(
       }
       else if (oper_type == Operator::ANY_TYPE)
       {
-         const Operator *out_R = fes_out.GetRestrictionOperator();
          t_oper.Reset(new TripleProductOperator(
                          out_R, &oper, fes_in.GetProlongationMatrix(),
                          false, false, false));

fem/fespace.hpp

@@ -330,18 +330,6 @@ public:
    virtual const Operator *GetProlongationMatrix() const
    { return GetConformingProlongation(); }
 
-   /// Return an operator that performs the transpose of GetRestrictionOperator
-   /** The returned operator is owned by the FiniteElementSpace. In serial this
-       is the same as GetProlongationMatrix() */
-   virtual const Operator *GetRestrictionTransposeOperator() const
-   { return GetConformingProlongation(); }
-
-   /// An abstract operator that performs the same action as GetRestrictionMatrix
-   /** In some cases this is an optimized matrix-free implementation. The
-       returned operator is owned by the FiniteElementSpace. */
-   virtual const Operator *GetRestrictionOperator() const
-   { return GetConformingRestriction(); }
-
    /// The returned SparseMatrix is owned by the FiniteElementSpace.
    virtual const SparseMatrix *GetRestrictionMatrix() const
    { return GetConformingRestriction(); }

fem/geom.cpp

@@ -1165,31 +1165,24 @@ RefinedGeometry * GeometryRefiner::Refine(Geometry::Type Geom,
          Array<int> vi((n+1)*(n+1)*(n+1));
          vi = -1;
          m = 0;
-
-         // 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++)
+         for (k = 0; k <= n; k++)
+            for (j = 0; j <= k; j++)
+               for (i = 0; i <= j; i++)
                {
                   IntegrationPoint &ip = RG->RefPts.IntPoint(m);
-                  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;
+                  // 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;
                   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");
@@ -1276,9 +1269,18 @@ RefinedGeometry * GeometryRefiner::Refine(Geometry::Type Geom,
                for (i = 0; i <= n-j; i++, l++)
                {
                   IntegrationPoint &ip = RG->RefPts.IntPoint(l);
-                  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];
+                  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];
+                  }
                   m++;
                }
          if (m != (n+1)*(n+1)*(n+2)/2)

fem/lininteg.hpp

@@ -406,10 +406,6 @@ 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,92 +128,6 @@ 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,80 +130,6 @@ 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.cpp

@@ -565,38 +565,6 @@ HypreParMatrix* ParDiscreteLinearOperator::ParallelAssemble() const
    return RAP;
 }
 
-void ParDiscreteLinearOperator::ParallelAssemble(OperatorHandle &A)
-{
-   // construct the rectangular block-diagonal matrix dA
-   OperatorHandle dA(A.Type());
-   dA.MakeRectangularBlockDiag(domain_fes->GetComm(),
-                               range_fes->GlobalVSize(),
-                               domain_fes->GlobalVSize(),
-                               range_fes->GetDofOffsets(),
-                               domain_fes->GetDofOffsets(),
-                               mat);
-
-   OperatorHandle R_test_transpose(A.Type()), P_trial(A.Type());
-
-   // TODO - construct the Dof_TrueDof_Matrix directly in the required format.
-   R_test_transpose.ConvertFrom(range_fes->Dof_TrueDof_Matrix());
-   P_trial.ConvertFrom(domain_fes->Dof_TrueDof_Matrix());
-
-   A.MakeRAP(R_test_transpose, dA, P_trial);
-}
-
-void ParDiscreteLinearOperator::FormRectangularSystemMatrix(OperatorHandle &A)
-{
-   if (ext)
-   {
-      Array<int> empty;
-      ext->FormRectangularSystemOperator(empty, empty, A);
-      return;
-   }
-
-   mfem_error("not implemented!");
-}
-
 void ParDiscreteLinearOperator::GetParBlocks(Array2D<HypreParMatrix *> &blocks)
 const
 {

fem/pbilinearform.hpp

@@ -160,9 +160,6 @@ public:
    /// Get the parallel finite element space prolongation matrix
    virtual const Operator *GetProlongation() const
    { return pfes->GetProlongationMatrix(); }
-   /// Get the transpose of GetRestriction, useful for matrix-free RAP
-   virtual const Operator *GetRestrictionTranspose() const
-   { return pfes->GetRestrictionTransposeOperator(); }
    /// Get the parallel finite element space restriction matrix
    virtual const Operator *GetRestriction() const
    { return pfes->GetRestrictionMatrix(); }
@@ -249,9 +246,6 @@ 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
@@ -307,20 +301,10 @@ public:
    /// Returns the matrix "assembled" on the true dofs
    HypreParMatrix *ParallelAssemble() const;
 
-   /** @brief Returns the matrix assembled on the true dofs, i.e.
-       @a A = R_test A_local P_trial, in the format (type id) specified by
-       @a A. */
-   void ParallelAssemble(OperatorHandle &A);
-
    /** Extract the parallel blocks corresponding to the vector dimensions of the
        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);
-
    virtual ~ParDiscreteLinearOperator() { }
 };
 

fem/pfespace.cpp

@@ -101,8 +101,6 @@ void ParFiniteElementSpace::ParInit(ParMesh *pm)
 
    P = NULL;
    Pconf = NULL;
-   Rconf = NULL;
-   R_transpose = NULL;
    R = NULL;
 
    num_face_nbr_dofs = -1;
@@ -929,45 +927,6 @@ const Operator *ParFiniteElementSpace::GetProlongationMatrix() const
    }
 }
 
-const Operator *ParFiniteElementSpace::GetRestrictionOperator() const
-{
-   if (Conforming())
-   {
-      if (Rconf) { return Rconf; }
-
-      if (NRanks == 1)
-      {
-         R_transpose = new IdentityOperator(GetTrueVSize());
-      }
-      else
-      {
-         if (!Device::Allows(Backend::DEVICE_MASK))
-         {
-            R_transpose = new ConformingProlongationOperator(*this, true);
-         }
-         else
-         {
-            R_transpose =
-               new DeviceConformingProlongationOperator(*this, true);
-         }
-      }
-      Rconf = new TransposeOperator(R_transpose);
-      return Rconf;
-   }
-   else
-   {
-      Dof_TrueDof_Matrix();
-      R_transpose = new TransposeOperator(R);
-      return R;
-   }
-}
-
-const Operator *ParFiniteElementSpace::GetRestrictionTransposeOperator() const
-{
-   GetRestrictionOperator();
-   return R_transpose;
-}
-
 void ParFiniteElementSpace::ExchangeFaceNbrData()
 {
    if (num_face_nbr_dofs >= 0) { return; }
@@ -2881,8 +2840,6 @@ void ParFiniteElementSpace::Destroy()
 
    delete P; P = NULL;
    delete Pconf; Pconf = NULL;
-   delete Rconf; Rconf = NULL;
-   delete R_transpose; R_transpose = NULL;
    delete R; R = NULL;
 
    delete gcomm; gcomm = NULL;
@@ -3008,12 +2965,12 @@ void ParFiniteElementSpace::Update(bool want_transform)
    }
 }
 
+
 ConformingProlongationOperator::ConformingProlongationOperator(
-   const ParFiniteElementSpace &pfes, bool local_)
+   const ParFiniteElementSpace &pfes)
    : Operator(pfes.GetVSize(), pfes.GetTrueVSize()),
      external_ldofs(),
-     gc(pfes.GroupComm()),
-     local(local_)
+     gc(pfes.GroupComm())
 {
    MFEM_VERIFY(pfes.Conforming(), "");
    const Table &group_ldof = gc.GroupLDofTable();
@@ -3062,14 +3019,7 @@ void ConformingProlongationOperator::Mult(const Vector &x, Vector &y) const
    const int m = external_ldofs.Size();
 
    const int in_layout = 2; // 2 - input is ltdofs array
-   if (local)
-   {
-      y = 0.0;
-   }
-   else
-   {
-      gc.BcastBegin(const_cast<double*>(xdata), in_layout);
-   }
+   gc.BcastBegin(const_cast<double*>(xdata), in_layout);
 
    int j = 0;
    for (int i = 0; i < m; i++)
@@ -3081,10 +3031,7 @@ void ConformingProlongationOperator::Mult(const Vector &x, Vector &y) const
    std::copy(xdata+j-m, xdata+Width(), ydata+j);
 
    const int out_layout = 0; // 0 - output is ldofs array
-   if (!local)
-   {
-      gc.BcastEnd(ydata, out_layout);
-   }
+   gc.BcastEnd(ydata, out_layout);
 }
 
 void ConformingProlongationOperator::MultTranspose(
@@ -3097,10 +3044,7 @@ void ConformingProlongationOperator::MultTranspose(
    double *ydata = y.HostWrite();
    const int m = external_ldofs.Size();
 
-   if (!local)
-   {
-      gc.ReduceBegin(xdata);
-   }
+   gc.ReduceBegin(xdata);
 
    int j = 0;
    for (int i = 0; i < m; i++)
@@ -3112,18 +3056,13 @@ void ConformingProlongationOperator::MultTranspose(
    std::copy(xdata+j, xdata+Height(), ydata+j-m);
 
    const int out_layout = 2; // 2 - output is an array on all ltdofs
-   if (!local)
-   {
-      gc.ReduceEnd<double>(ydata, out_layout, GroupCommunicator::Sum);
-   }
+   gc.ReduceEnd<double>(ydata, out_layout, GroupCommunicator::Sum);
 }
 
 DeviceConformingProlongationOperator::DeviceConformingProlongationOperator(
-   const ParFiniteElementSpace &pfes,
-   bool local_) :
+   const ParFiniteElementSpace &pfes) :
    ConformingProlongationOperator(pfes),
-   mpi_gpu_aware(Device::GetGPUAwareMPI()),
-   local(local_)
+   mpi_gpu_aware(Device::GetGPUAwareMPI())
 {
    MFEM_ASSERT(pfes.Conforming(), "internal error");
    const SparseMatrix *R = pfes.GetRestrictionMatrix();
@@ -3240,42 +3179,32 @@ void DeviceConformingProlongationOperator::Mult(const Vector &x,
                                                 Vector &y) const
 {
    const GroupTopology &gtopo = gc.GetGroupTopology();
+   BcastBeginCopy(x); // copy to 'shr_buf'
    int req_counter = 0;
-   if (local)
+   for (int nbr = 1; nbr < gtopo.GetNumNeighbors(); nbr++)
    {
-      y = 0.0;
-   }
-   else
-   {
-      BcastBeginCopy(x); // copy to 'shr_buf'
-      for (int nbr = 1; nbr < gtopo.GetNumNeighbors(); nbr++)
+      const int send_offset = shr_buf_offsets[nbr];
+      const int send_size = shr_buf_offsets[nbr+1] - send_offset;
+      if (send_size > 0)
       {
-         const int send_offset = shr_buf_offsets[nbr];
-         const int send_size = shr_buf_offsets[nbr+1] - send_offset;
-         if (send_size > 0)
-         {
-            auto send_buf = mpi_gpu_aware ? shr_buf.Read() : shr_buf.HostRead();
-            MPI_Isend(send_buf + send_offset, send_size, MPI_DOUBLE,
-                      gtopo.GetNeighborRank(nbr), 41822,
-                      gtopo.GetComm(), &requests[req_counter++]);
-         }
-         const int recv_offset = ext_buf_offsets[nbr];
-         const int recv_size = ext_buf_offsets[nbr+1] - recv_offset;
-         if (recv_size > 0)
-         {
-            auto recv_buf = mpi_gpu_aware ? ext_buf.Write() : ext_buf.HostWrite();
-            MPI_Irecv(recv_buf + recv_offset, recv_size, MPI_DOUBLE,
-                      gtopo.GetNeighborRank(nbr), 41822,
-                      gtopo.GetComm(), &requests[req_counter++]);
-         }
+         auto send_buf = mpi_gpu_aware ? shr_buf.Read() : shr_buf.HostRead();
+         MPI_Isend(send_buf + send_offset, send_size, MPI_DOUBLE,
+                   gtopo.GetNeighborRank(nbr), 41822,
+                   gtopo.GetComm(), &requests[req_counter++]);
+      }
+      const int recv_offset = ext_buf_offsets[nbr];
+      const int recv_size = ext_buf_offsets[nbr+1] - recv_offset;
+      if (recv_size > 0)
+      {
+         auto recv_buf = mpi_gpu_aware ? ext_buf.Write() : ext_buf.HostWrite();
+         MPI_Irecv(recv_buf + recv_offset, recv_size, MPI_DOUBLE,
+                   gtopo.GetNeighborRank(nbr), 41822,
+                   gtopo.GetComm(), &requests[req_counter++]);
       }
    }
    BcastLocalCopy(x, y);
-   if (!local)
-   {
-      MPI_Waitall(req_counter, requests, MPI_STATUSES_IGNORE);
-      BcastEndCopy(y); // copy from 'ext_buf'
-   }
+   MPI_Waitall(req_counter, requests, MPI_STATUSES_IGNORE);
+   BcastEndCopy(y); // copy from 'ext_buf'
 }
 
 DeviceConformingProlongationOperator::~DeviceConformingProlongationOperator()
@@ -3338,38 +3267,32 @@ void DeviceConformingProlongationOperator::MultTranspose(const Vector &x,
                                                          Vector &y) const
 {
    const GroupTopology &gtopo = gc.GetGroupTopology();
+   ReduceBeginCopy(x); // copy to 'ext_buf'
    int req_counter = 0;
-   if (!local)
+   for (int nbr = 1; nbr < gtopo.GetNumNeighbors(); nbr++)
    {
-      ReduceBeginCopy(x); // copy to 'ext_buf'
-      for (int nbr = 1; nbr < gtopo.GetNumNeighbors(); nbr++)
+      const int send_offset = ext_buf_offsets[nbr];
+      const int send_size = ext_buf_offsets[nbr+1] - send_offset;
+      if (send_size > 0)
       {
-         const int send_offset = ext_buf_offsets[nbr];
-         const int send_size = ext_buf_offsets[nbr+1] - send_offset;
-         if (send_size > 0)
-         {
-            auto send_buf = mpi_gpu_aware ? ext_buf.Read() : ext_buf.HostRead();
-            MPI_Isend(send_buf + send_offset, send_size, MPI_DOUBLE,
-                      gtopo.GetNeighborRank(nbr), 41823,
-                      gtopo.GetComm(), &requests[req_counter++]);
-         }
-         const int recv_offset = shr_buf_offsets[nbr];
-         const int recv_size = shr_buf_offsets[nbr+1] - recv_offset;
-         if (recv_size > 0)
-         {
-            auto recv_buf = mpi_gpu_aware ? shr_buf.Write() : shr_buf.HostWrite();
-            MPI_Irecv(recv_buf + recv_offset, recv_size, MPI_DOUBLE,
-                      gtopo.GetNeighborRank(nbr), 41823,
-                      gtopo.GetComm(), &requests[req_counter++]);
-         }
+         auto send_buf = mpi_gpu_aware ? ext_buf.Read() : ext_buf.HostRead();
+         MPI_Isend(send_buf + send_offset, send_size, MPI_DOUBLE,
+                   gtopo.GetNeighborRank(nbr), 41823,
+                   gtopo.GetComm(), &requests[req_counter++]);
+      }
+      const int recv_offset = shr_buf_offsets[nbr];
+      const int recv_size = shr_buf_offsets[nbr+1] - recv_offset;
+      if (recv_size > 0)
+      {
+         auto recv_buf = mpi_gpu_aware ? shr_buf.Write() : shr_buf.HostWrite();
+         MPI_Irecv(recv_buf + recv_offset, recv_size, MPI_DOUBLE,
+                   gtopo.GetNeighborRank(nbr), 41823,
+                   gtopo.GetComm(), &requests[req_counter++]);
       }
    }
    ReduceLocalCopy(x, y);
-   if (!local)
-   {
-      MPI_Waitall(req_counter, requests, MPI_STATUSES_IGNORE);
-      ReduceEndAssemble(y); // assemble from 'shr_buf'
-   }
+   MPI_Waitall(req_counter, requests, MPI_STATUSES_IGNORE);
+   ReduceEndAssemble(y); // assemble from 'shr_buf'
 }
 
 } // namespace mfem

fem/pfespace.hpp

@@ -75,12 +75,6 @@ private:
 
    /// The (block-diagonal) matrix R (restriction of dof to true dof). Owned.
    mutable SparseMatrix *R;
-   /// Optimized action-only restriction operator for conforming meshes. Owned.
-   mutable Operator *Rconf;
-   /** Transpose of R or Rconf. For conforming mesh, this is a matrix-free
-       (Device)ConformingProlongationOperator, for a non-conforming mesh
-       this is a TransposeOperator wrapping R. */
-   mutable Operator *R_transpose;
 
    ParNURBSExtension *pNURBSext() const
    { return dynamic_cast<ParNURBSExtension *>(NURBSext); }
@@ -347,16 +341,6 @@ public:
    HYPRE_Int GetMyTDofOffset() const;
 
    virtual const Operator *GetProlongationMatrix() const;
-   /** @brief Return logical transpose of restriction matrix, but in
-       non-assembled optimized matrix-free form.
-
-       The implementation is like GetProlongationMatrix, but it sets local
-       DOFs to the true DOF values if owned locally, otherwise zero. */
-   virtual const Operator *GetRestrictionTransposeOperator() const;
-   /** Get an Operator that performs the action of GetRestrictionMatrix(),
-       but potentially with a non-assembled optimized matrix-free
-       implementation. */
-   virtual const Operator *GetRestrictionOperator() const;
    /// Get the R matrix which restricts a local dof vector to true dof vector.
    virtual const SparseMatrix *GetRestrictionMatrix() const
    { Dof_TrueDof_Matrix(); return R; }
@@ -411,11 +395,9 @@ class ConformingProlongationOperator : public Operator
 protected:
    Array<int> external_ldofs;
    const GroupCommunicator &gc;
-   bool local;
 
 public:
-   ConformingProlongationOperator(const ParFiniteElementSpace &pfes,
-                                  bool local_=false);
+   ConformingProlongationOperator(const ParFiniteElementSpace &pfes);
 
    virtual void Mult(const Vector &x, Vector &y) const;
 
@@ -434,8 +416,6 @@ protected:
    Array<int> ltdof_ldof, unq_ltdof;
    Array<int> unq_shr_i, unq_shr_j;
    MPI_Request *requests;
-   bool local;
-
    // Kernel: copy ltdofs from 'src' to 'shr_buf' - prepare for send.
    //         shr_buf[i] = src[shr_ltdof[i]]
    void BcastBeginCopy(const Vector &src) const;
@@ -461,8 +441,7 @@ protected:
    void ReduceEndAssemble(Vector &dst) const;
 
 public:
-   DeviceConformingProlongationOperator(const ParFiniteElementSpace &pfes,
-                                        bool local_=false);
+   DeviceConformingProlongationOperator(const ParFiniteElementSpace &pfes);
 
    virtual ~DeviceConformingProlongationOperator();
 

fem/pprmnonlinearform.cpp

@@ -1,363 +0,0 @@
-// 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

@@ -1,105 +0,0 @@
-// 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

@@ -1,231 +0,0 @@
-// 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
-
-

fem/restriction.cpp

@@ -195,31 +195,6 @@ void ElementRestriction::MultTransposeUnsigned(const Vector& x, Vector& y) const
    });
 }
 
-void ElementRestriction::MultLeftInverse(const Vector& x, Vector& y) const
-{
-   // Assumes all elements have the same number of dofs
-   const int nd = dof;
-   const int vd = vdim;
-   const bool t = byvdim;
-   auto d_offsets = offsets.Read();
-   auto d_indices = indices.Read();
-   auto d_x = Reshape(x.Read(), nd, vd, ne);
-   auto d_y = Reshape(y.Write(), t?vd:ndofs, t?ndofs:vd);
-   MFEM_FORALL(i, ndofs,
-   {
-      const int nextOffset = d_offsets[i + 1];
-      for (int c = 0; c < vd; ++c)
-      {
-         double dofValue = 0;
-         const int j = nextOffset - 1;
-         const int idx_j = (d_indices[j] >= 0) ? d_indices[j] : -1 - d_indices[j];
-         dofValue = (d_indices[j] >= 0) ? d_x(idx_j % nd, c, idx_j / nd) :
-         -d_x(idx_j % nd, c, idx_j / nd);
-         d_y(t?c:i,t?i:c) = dofValue;
-      }
-   });
-}
-
 void ElementRestriction::BooleanMask(Vector& y) const
 {
    // Assumes all elements have the same number of dofs

fem/restriction.hpp

@@ -57,10 +57,6 @@ public:
    /// Compute MultTranspose without applying signs based on DOF orientations.
    void MultTransposeUnsigned(const Vector &x, Vector &y) const;
 
-   /// Compute MultTranspose by setting (rather than adding) element
-   /// contributions; this is a left inverse of the Mult() operation
-   void MultLeftInverse(const Vector &x, Vector &y) const;
-
    /// @brief Fills the E-vector y with `boolean` values 0.0 and 1.0 such that each
    /// each entry of the L-vector is uniquely represented in `y`.
    /** This means, the sum of the E-vector `y` is equal to the sum of the

general/CMakeLists.txt

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

linalg/hypre.cpp

@@ -128,15 +128,6 @@ 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);
@@ -935,15 +926,6 @@ 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
@@ -984,46 +966,6 @@ 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)
 {
@@ -1627,36 +1569,6 @@ 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,
@@ -2905,11 +2817,6 @@ 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);
@@ -3713,116 +3620,6 @@ 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,14 +81,6 @@ 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
@@ -121,9 +113,6 @@ 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. */
@@ -239,24 +228,15 @@ 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();
-      WrapHypreParCSRMatrix(a, owner);
+      A = a;
+      if (!owner) { ParCSROwner = 0; }
+      height = GetNumRows();
+      width = GetNumCols();
    }
 
    /// Creates block-diagonal square parallel matrix.
@@ -407,13 +387,6 @@ 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. */
@@ -424,13 +397,6 @@ 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
    {
@@ -583,23 +549,6 @@ 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. */
@@ -700,12 +649,11 @@ 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, OPFS = 10, Chebyshev = 16, Taubin = 1001, FIR = 1002
+               GS = 6, Chebyshev = 16, Taubin = 1001, FIR = 1002
              };
 
    HypreSmoother();
@@ -819,28 +767,6 @@ 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
 {
@@ -918,7 +844,6 @@ 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);
@@ -1203,94 +1128,9 @@ 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; }
 
@@ -1299,8 +1139,6 @@ public:
    virtual HYPRE_PtrToParSolverFcn SolveFcn() const
    { return (HYPRE_PtrToParSolverFcn) HYPRE_BoomerAMGSolve; }
 
-   using HypreSolver::Mult;
-
    virtual ~HypreBoomerAMG();
 };
 

linalg/operator.hpp

@@ -37,8 +37,7 @@ protected:
       const Array<int> &test_tdof_list,
       RectangularConstrainedOperator* &Aout);
 
-   /** @brief Returns RAP Operator of this, using input/output Prolongation matrices
-       @a Pi corresponds to "P", @a Po corresponds to "Rt" */
+   /// Returns RAP Operator of this, taking in input/output Prolongation matrices
    Operator *SetupRAP(const Operator *Pi, const Operator *Po);
 
 public:
@@ -113,11 +112,6 @@ public:
    {
       return GetProlongation(); // Assume square unless specialized
    }
-   /** @brief Transpose of GetOutputRestriction, directly available in this
-       form to facilitate matrix-free RAP-type operators.
-
-       `NULL` means identity. */
-   virtual const Operator *GetOutputRestrictionTranspose() const { return NULL; }
    /** @brief Restriction operator from output vectors for the operator to linear
        algebra (linear system) vectors. `NULL` means identity. */
    virtual const Operator *GetOutputRestriction() const
@@ -657,7 +651,6 @@ public:
    virtual void SetOperator(const Operator &op) = 0;
 };
 
-
 /// Identity Operator I: x -> x.
 class IdentityOperator : public Operator
 {

linalg/solvers.cpp

@@ -1733,6 +1733,323 @@ void NewtonSolver::AdaptiveLinRtolPostSolve(const Vector &x,
    }
 }
 
+void AndersonAcceleration::QRdelete(std::deque<Vector *> &Q, DenseMatrix &R) const
+{
+   Vector temp(Q[0]->Size());
+   for (int i=0; i<(maxVecs-1); i++) {
+      double d = sqrt( R(i,i+1)*R(i,i+1) + R(i+1,i+1)*R(i+1,i+1) );
+      double c = R(i,i+1) / d;
+      double s = R(i+1,i+1) / d;
+      R(i,i+1) = d;
+      R(i+1,i+1) = 0;
+
+      if (i < (maxVecs-2)) {
+         for (int j=(i+2); j<maxVecs; j++) {
+            d = c*R(i,j) + s*R(i+1,j);
+            R(i+1,j) = -s*R(i,j) + c*R(i+1,j);
+            R(i,j) = d;
+         }
+      }
+      // temp = c*Q[i] + s*Q[i+1];
+      add(c, *(Q[i]), s, *(Q[i+1]), temp);
+      // Q[i+1] = -s*Q[i] + c*Q[i+1];
+      *(Q[i+1]) *= c;
+      Q[i+1] -> Add(-s, *(Q[i]));
+      *(Q[i]) = temp;
+   }
+   
+   // Shift Q <- Q[:,0:(m-2)], i.e., delete last column of Q
+   delete Q.back();
+   Q.pop_back();
+
+   // Shift columns of R to the left by one, R = R[0:(m−2), 1:(m-1)]
+   for (int j=1; j<maxVecs; j++) {
+      for (int i=0; i<maxVecs; i++) {
+         R(i,j-1) = R(i,j);
+      }
+   }
+}
+
+void AndersonAcceleration::SetOperator(const Operator &op)
+{
+   oper = &op;
+   height = op.Height();
+   width = op.Width();
+   MFEM_ASSERT(height == width, "square Operator is required.");
+}
+
+void AndersonAcceleration::FixedPointMult(const Vector &b,
+   const Vector &x, Vector &y, Vector &r) const
+{
+   // Assume Operator represents a fixed-point operator *with
+   // right-hand side*, so we iterate x_{k+1} = M^{-1}G(x_k)
+   if (isFixedPointOp) {
+      oper->Mult(x, y);
+      if (prec) {
+         prec->Mult(y,r);
+         y = r;
+         r -= x;
+      }
+      else {
+         r = y;
+         r -= x;
+      }
+   }
+   // Otherwise add x to y and not to residual r = y - x
+   else {
+      oper->Mult(x, r);
+      r *= -1;
+      r += b;
+      if (prec) {
+         prec->Mult(r,y);
+         r = y;
+         y += x;
+      }
+      else {
+         y = r;
+         y += x;         
+      }
+   }
+}
+
+void AndersonAcceleration::Mult(const Vector &b, Vector &x) const
+{
+   MFEM_ASSERT(oper != NULL, "the Operator is not set (use SetOperator).");
+  
+   int n = width;
+   int numVecs = 0;
+   double resid, norm_df;
+   double min_diag = 1e-13;
+   final_norm = -1;
+
+   // Check vector is initialized, set to zero for
+   // iterative_mode = false
+   if (x.Size() != n) {
+      x.SetSize(n);
+      x = 0.0;
+   }
+   else if (!iterative_mode) {
+      x = 0.0;
+   }
+
+   // Storage containers for acceleration
+   std::deque<Vector *> G;
+   std::deque<Vector *> Q;
+   DenseMatrix R(maxVecs);
+   R = 0.0;
+   Vector g_old(n);
+   Vector g_current(n);
+   Vector f_old(n);
+   Vector f_current(n);
+   Vector gamma(maxVecs);
+   Vector rhs(maxVecs);
+   Vector correction;
+   if (omega > 0 && std::abs(omega -  1) > 1e-14) {
+      correction.SetSize(n);
+   }
+   Vector *dg;
+   Vector *df;
+
+   // Loop over AA iterations
+   int k;
+   for (k=0; k<max_iter; k++) {
+
+      // Compute g_current = G(x), f_current = G(x) - x
+      this->FixedPointMult(b, x, g_current, f_current); 
+
+      // Check norm of current approximation to fixed point G(u) = u
+      resid = Norm(f_current);
+      MFEM_ASSERT(IsFinite(resid), "||G(u) - u|| = " << resid);
+      if (print_level == 1)
+      {
+         mfem::out << "  Iteration : " << setw(3) << k
+                 << "  ||G(u) - u|| = " << resid << endl;
+      }
+      
+      // Set stopping tolerance on first iteration.
+      if (final_norm < 0) {
+         final_norm = std::max(rel_tol*resid, abs_tol);
+      }
+
+      // Check for convergence
+      if (resid <= final_norm)
+      {
+         final_norm = resid;
+         final_iter = k;
+         converged = 1;
+         goto finish;
+      }
+
+      // Start Anderson Acceleration after AAstart FP iterations
+      if (k > AAstart) {
+         // df = f_current - f_old; 
+         df = new Vector(n);
+         add(1.0, f_current, -1.0, f_old, *df);
+         // dg = g_current - g_old;
+         dg = new Vector(n);
+         add(1.0, g_current, -1.0, g_old, *dg);
+      
+         if (numVecs < maxVecs) {
+            G.push_back(dg);
+         }
+         else {
+            delete G[0];
+            G.pop_front();
+            G.push_back(dg);
+         }
+         numVecs++;
+         dg = NULL;
+      }
+      
+      f_old = f_current;
+      g_old = g_current;
+      
+      // First iteration or initial fixed-point iterations
+      if (numVecs == 0) {
+         x = g_current;
+         continue;
+      }
+
+      // All later iterations: orthogonalize and find best approximation
+      if (numVecs == 1) {
+         norm_df = Norm(*df);
+         MFEM_ASSERT(IsFinite(norm_df), "norm_df = " << norm_df);
+         (*df) /= norm_df;
+         Q.push_back(df);
+         R(0,0) = norm_df;
+         df = NULL;
+      }
+      else {
+         // Remove first column in basis F and R, reorthogonalize
+         if (numVecs > maxVecs) {
+            this->QRdelete(Q, R);
+            numVecs--;
+         }
+         // Compute last column of R
+         for (int i=0; i<(numVecs-1); i++) {
+            R(i,numVecs-1) = Dot(*(Q[i]), *df);
+            // df -= R(i,numVecs-1) * Q[i]
+            df -> Add(-R(i,numVecs-1), *(Q[i]));
+         }
+         norm_df = Norm(*df);
+         MFEM_ASSERT(IsFinite(norm_df), "norm_df = " << norm_df);
+         (*df) /= norm_df;
+         Q.push_back(df);
+         R(numVecs-1, numVecs-1) = norm_df;
+         df = NULL;
+      }
+
+      // Back solve for new weights, R\gamma = Q^T * f_current
+      rhs = 0.0;
+      gamma = 0.0;
+      for (int i=0; i<numVecs; i++) {
+         rhs(i) = Dot(*(Q[i]), f_current);   // Form right hand side
+      }
+      for (int i=(numVecs-1); i>=0; i--) {
+         double temp = rhs(i);
+         for (int j=(i+1); j<numVecs; j++) {
+            temp -= R(i,j)*gamma(j);
+         }
+         if (std::abs(R(i,i)) < min_diag) {
+            gamma(i) = 0.0;
+            std::cout << "Diagonal of R -- " << R(i,i) << " ~ 0.\n";
+         }
+         else {
+            gamma(i) = temp / R(i,i);            
+         }
+      }
+
+      /// DEBUG --> test backsolve
+      Vector test(numVecs);
+      for (int i=0; i<numVecs; i++) {
+         test(i) = 0;
+         for (int j=i; j<numVecs; j++) {
+            test(i) += R(i,j) * gamma(j);
+         }
+         if (std::abs(rhs(i) - test(i)) > 1e-10) {
+            std::cout << "Bad solve! Err = " << rhs(i) - test(i) << "\n";
+         }
+      }
+
+      // Compute updated solution x = g_current − G*\gamma 
+      x = g_current;
+      for (int i=0; i<numVecs; i++) {
+         // x -= gamma(i)*G[i]
+         x.Add(-gamma(i), *(G[i]));
+      }
+
+      // Apply damped iteration for \omega \in (0,1),
+      //   x -= (1−omega) * (f_current − Q*R*gamma);
+      if (omega > 0 && std::abs(omega -  1) > 1e-14) {
+         // Redefine rhs = R*gamma
+         for(int i=0; i<numVecs; i++) {
+            rhs(i) = 0;
+            for (int j=i; j<numVecs; j++) {
+               rhs(i) += R(i,j)*gamma(j);
+            }
+         }
+         correction = f_current;
+         for (int i=0; i<numVecs; i++) {
+            // correction -= rhs(i)*Q[i]
+            correction.Add(-rhs(i), *(Q[i]));
+         }
+         // x -= (1 - omega) * correction;
+         x.Add( -(1 - omega), correction);
+      }
+
+      // Restart AA minimization by eliminating all vectors but the most recent
+      if (restart && (numVecs == maxVecs)) {
+         for (int i=0; i<(maxVecs-1); i++) {
+            delete G[0];
+            G.pop_front();          
+            delete Q[0];
+            Q.pop_front();
+         }
+         R = 0.0;
+         R(0,0) = norm_df;
+         numVecs = 1;
+         if (print_level == 1)
+         {
+            mfem::out << "Restarting..." << '\n';
+         }
+      }
+   }
+
+   // Compute final residual, save counts for solve
+   this->FixedPointMult(b, x, g_current, f_current); 
+   resid = Norm(f_current);
+   MFEM_ASSERT(IsFinite(resid), "||G(u) - u|| = " << resid);
+   final_norm = resid;
+   final_iter = max_iter;
+   if (resid <= final_norm) converged = 1;
+   else converged = 0;
+
+finish:
+   if (print_level == 3)
+   {
+      mfem::out << "  Iteration : " << setw(3) << k
+              << "  ||G(u) - u|| = " << resid << endl;
+   }
+   else if (print_level == 2)
+   {
+      mfem::out << "Anderson Acceleration: Number of iterations: " << final_iter << '\n';
+   }
+   if (print_level >= 0 && !converged)
+   {
+      mfem::out << "Anderson Acceleration: No convergence!\n";
+   }
+
+   // Cleanup pointers
+   for (int i=0; i<numVecs; i++) {
+      delete G[0];
+      G.pop_front();          
+      delete Q[0];
+      Q.pop_front();
+   }
+   delete dg;
+   delete df;
+}
+
 void LBFGSSolver::Mult(const Vector &b, Vector &x) const
 {
    MFEM_VERIFY(oper != NULL, "the Operator is not set (use SetOperator).");

linalg/solvers.hpp

@@ -14,6 +14,7 @@
 
 #include "../config/config.hpp"
 #include "densemat.hpp"
+#include <deque>
 
 #ifdef MFEM_USE_MPI
 #include <mpi.h>
@@ -489,6 +490,64 @@ public:
                            const double gamma = 1.0);
 };
 
+
+/// Nonlinear Anderson Acceleration
+class AndersonAcceleration : public IterativeSolver
+{
+protected:
+   int maxVecs; // see SetKDim()
+   int AAstart;
+   bool restart;
+   bool isFixedPointOp;
+   double omega;
+
+   /// Apply fixed-point Mult and compute residual.
+   //  For !isFixedPointOp:
+   //    y <-- G(x) = x + A(x) - b and r <-- G(x) - x = A(x) - b
+   //  For isFixedPointOp:
+   //    y <-- G(x) = A(x) - b and r <-- G(x) - x = A(x) - b - x
+   void FixedPointMult(const Vector &b, const Vector &x, Vector &y, Vector &r) const;
+
+   // Helper function for Anderson Acceleration
+   void QRdelete(std::deque<Vector *> &Q, DenseMatrix &R) const;
+
+public:
+   AndersonAcceleration() : maxVecs(25), AAstart(0), omega(1),
+      restart(false), isFixedPointOp(false) { }
+
+#ifdef MFEM_USE_MPI
+   AndersonAcceleration(MPI_Comm _comm) : IterativeSolver(_comm),
+      maxVecs(25), AAstart(0), omega(1), restart(false),
+      isFixedPointOp(false) { }
+#endif
+
+   /// Boolean describing whether the action of the operator is such
+   //  that we want to solve G(x) = x (true) or G(x) = 0 (false) for
+   //  zero right-hand side vector b passed into Mult(). Default is
+   //  false in construction of class.
+   void IsFixedPointOperator(bool isFixedPointOp_)
+   { isFixedPointOp = isFixedPointOp_; }
+
+   /// Maximum number of vectors to store in Krylov-like space
+   void SetKDim(int dim) { maxVecs = dim; }
+
+   /// Number of fixed-point iterations to do before starting AA
+   void SetAAStart(int start_) { AAstart = start_; }
+
+   /// Boolean to restart, that is, erase entire space after maxVecs
+   //  are stored (AAstart=true) or use a sliding space (AAstart=false)
+   //  where one vector is deleted to make room for a new one.
+   void SetRestart(bool restart_) { restart = restart_; }
+
+   /// Set relaxation weight
+   void SetWeight(double omega_) { omega = omega_; }
+
+   virtual void SetOperator(const Operator &op);
+
+   virtual void Mult(const Vector &b, Vector &x) const;
+};
+
+
 /** L-BFGS method for solving F(x)=b for a given operator F, by minimizing
     the norm of F(x) - b. Requires only the action of the operator F. */
 class LBFGSSolver : public NewtonSolver

makefile

@@ -684,10 +684,9 @@ 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 := $(filter-out general/tinyxml2.cpp,$(wildcard $(FORMAT_FILES)))
+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")
 
 COUT_CERR_FILES = $(foreach dir,$(DIRS),$(dir)/*.[ch]pp)
 COUT_CERR_EXCLUDE = '^general/error\.cpp' '^general/globals\.[ch]pp'

mesh/mesh.cpp

@@ -3448,7 +3448,6 @@ 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);
@@ -3490,20 +3489,11 @@ void Mesh::Loader(std::istream &input, int generate_edges,
    {
       ReadTrueGridMesh(input);
    }
-   else if (mesh_type.rfind("# vtk DataFile Version") == 0)
+   else if (mesh_type == "# vtk DataFile Version 3.0" ||
+            mesh_type == "# vtk DataFile Version 2.0") // VTK
    {
-      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);
@@ -8948,7 +8938,7 @@ void Mesh::PrintVTK(std::ostream &out)
          const int nv = elements[i]->GetNVertices();
          out << nv;
          Geometry::Type geom = elements[i]->GetGeometryType();
-         const int *perm = VTKGeometry::VertexPermutation[geom];
+         const int *perm = (geom == Geometry::PRISM) ? vtk_prism_perm : NULL;
          for (int j = 0; j < nv; j++)
          {
             out << ' ' << v[perm ? perm[j] : j];
@@ -9033,9 +9023,35 @@ void Mesh::PrintVTK(std::ostream &out)
    for (int i = 0; i < NumOfElements; i++)
    {
       int vtk_cell_type = 5;
-      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]; }
+      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;
+         }
+      }
+
       out << vtk_cell_type << '\n';
    }
 
@@ -9253,7 +9269,7 @@ void Mesh::PrintVTU(std::ostream &out, int ref, VTKFormat format,
          {
             coff = coff+nv;
             offset.push_back(coff);
-            const int *p = VTKGeometry::VertexPermutation[geom];
+            const int *p = (geom == Geometry::PRISM) ? vtk_prism_perm : NULL;
             for (int k = 0; k < nv; k++, j++)
             {
                WriteBinaryOrASCII(out, buf, np + RG[p ? p[j] : j], " ", format);
@@ -9284,14 +9300,39 @@ 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_cell_type = vtk_geom_map[geom];
+      // 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;
+      }
 
       if (high_order_output)
       {
@@ -9440,7 +9481,21 @@ 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 = VTKGeometry::Map[geom];
+      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;
+      }
 
       for (int j = 0; j < RG.Size(); j += nv)
       {

mesh/mesh.hpp

@@ -239,15 +239,8 @@ 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,13 +12,10 @@
 #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"
@@ -376,398 +373,6 @@ 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)
 {
@@ -777,6 +382,8 @@ 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);
@@ -805,78 +412,141 @@ 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 >> buff;
-      if (!input.good())
+      input >> np >> ws;
+      points.SetSize(3*np);
+      getline(input, buff); // "double"
+      for (i = 0; i < points.Size(); i++)
       {
-         MFEM_ABORT("VTK mesh does not have CELLS data!");
+         input >> points(i);
       }
    }
-   while (buff != "CELLS");
 
    // Read the cells
-   Array<int> cell_data, cell_offsets;
+   NumOfElements = n = 0;
+   Array<int> cells_data;
+   input >> ws >> buff;
    if (buff == "CELLS")
    {
-      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 >> NumOfElements >> n >> ws;
+      cells_data.SetSize(n);
+      for (i = 0; i < n; i++)
       {
-         int nv;
-         input >> nv;
-         cell_offsets[i] = offset + nv;
-         for (int j=0; j<nv; ++j)
-         {
-            input >> cell_data[offset + j];
-         }
-         offset += nv;
+         input >> cells_data[i];
       }
    }
 
    // Read the cell types
+   Dim = -1;
+   int order = -1;
    input >> ws >> buff;
-   Array<int> cell_types;
-   int ncells;
    if (buff == "CELL_TYPES")
    {
-      input >> ncells;
-      cell_types.Load(ncells, input);
+      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;
+      }
    }
 
-   // Read cell attributes
+   // Read 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 (buff.rfind("SCALARS material") == 0)
+      if (!strncmp(buff.c_str(), "SCALARS material", 16))
       {
          getline(input, buff); // "LOOKUP_TABLE default"
-         cell_attributes.Load(ncells, input);
+         for (i = 0; i < NumOfElements; i++)
+         {
+            input >> attr;
+            elements[i]->SetAttribute(attr);
+         }
       }
       else
       {
@@ -888,9 +558,138 @@ void Mesh::ReadVTKMesh(std::istream &input, int &curved, int &read_gf,
       input.seekg(sp);
    }
 
-   CreateVTKMesh(points, cell_data, cell_offsets, cell_types, cell_attributes,
-                 curved, read_gf, finalize_topo);
-} // end ReadVTKMesh
+   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;
+   }
+}
 
 void Mesh::ReadNURBSMesh(std::istream &input, int &curved, int &read_gf)
 {
@@ -916,8 +715,8 @@ void Mesh::ReadNURBSMesh(std::istream &input, int &curved, int &read_gf)
       NURBSext->SetCoordsFromPatches(*Nodes);
       own_nodes = 1;
       read_gf = 0;
-      spaceDim = Nodes->VectorDim();
-      for (int i = 0; i < spaceDim; i++)
+      int vd = Nodes->VectorDim();
+      for (int i = 0; i < vd; i++)
       {
          Vector vert_val;
          Nodes->GetNodalValues(vert_val, i+1);

mesh/vtk.cpp

@@ -18,137 +18,7 @@
 namespace mfem
 {
 
-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;
-}
+const int vtk_prism_perm[6] = {0, 2, 1, 3, 5, 4};
 
 int BarycentricToVTKTriangle(int *b, int ref)
 {
@@ -509,12 +379,15 @@ void CreateVTKElementConnectivity(Array<int> &con, Geometry::Type geom, int ref)
    {
       int idx = 0;
       int b[4];
-      for (b[2]=0; b[2]<=ref; b[2]++)
+      for (int k=0; k<=ref; k++)
       {
-         for (b[1]=0; b[1]<=ref-b[2]; b[1]++)
+         for (int j=0; j<=k; j++)
          {
-            for (b[0]=0; b[0]<=ref-b[1]-b[2]; b[0]++)
+            for (int i=0; i<=j; i++)
             {
+               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,46 +17,7 @@
 namespace mfem
 {
 
-// 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);
-};
+// Helpers for writing to the VTK format
 
 enum class VTKFormat
 {
@@ -65,6 +26,8 @@ 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,7 +42,6 @@ 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
@@ -51,17 +50,14 @@ 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
   fem/test_pa_coeff.cpp
   fem/test_pa_kernels.cpp
-  fem/test_pa_grad.cpp
-  fem/test_pa_idinterp.cpp
   fem/test_quadf_coef.cpp
   fem/test_quadraturefunc.cpp
-  fem/test_sum_bilin.cpp
+  fem/test_blocknonlinearform.cpp
   miniapps/test_sedov.cpp
 )
 

tests/unit/fem/test_lexicographic_ordering.cpp

@@ -1,64 +0,0 @@
-// 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_pa_grad.cpp

@@ -1,219 +0,0 @@
-// 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 "catch.hpp"
-#include "mfem.hpp"
-#include "unit_tests.hpp"
-
-using namespace mfem;
-
-double compare_pa_assembly(int dim, int num_elements, int order, bool transpose)
-{
-   Mesh * mesh;
-   if (num_elements == 0)
-   {
-      if (dim == 2)
-      {
-         mesh = new Mesh("../../data/star.mesh", order);
-      }
-      else
-      {
-         mesh = new Mesh("../../data/beam-hex.mesh", order);
-      }
-   }
-   else
-   {
-      if (dim == 2)
-      {
-         mesh = new Mesh(num_elements, num_elements, Element::QUADRILATERAL, true);
-      }
-      else
-      {
-         mesh = new Mesh(num_elements, num_elements, num_elements,
-                         Element::HEXAHEDRON, true);
-      }
-   }
-   FiniteElementCollection *h1_fec = new H1_FECollection(order, dim);
-   FiniteElementCollection *nd_fec = new ND_FECollection(order, dim);
-   FiniteElementSpace h1_fespace(mesh, h1_fec);
-   FiniteElementSpace nd_fespace(mesh, nd_fec);
-
-   DiscreteLinearOperator assembled_grad(&h1_fespace, &nd_fespace);
-   assembled_grad.AddDomainInterpolator(new GradientInterpolator);
-   const int skip_zeros = 1;
-   assembled_grad.Assemble(skip_zeros);
-   assembled_grad.Finalize(skip_zeros);
-   const SparseMatrix& assembled_grad_mat = assembled_grad.SpMat();
-
-   DiscreteLinearOperator pa_grad(&h1_fespace, &nd_fespace);
-   pa_grad.SetAssemblyLevel(AssemblyLevel::PARTIAL);
-   pa_grad.AddDomainInterpolator(new GradientInterpolator);
-   pa_grad.Assemble();
-   pa_grad.Finalize();
-
-   int insize, outsize;
-   if (transpose)
-   {
-      insize = nd_fespace.GetVSize();
-      outsize = h1_fespace.GetVSize();
-   }
-   else
-   {
-      insize = h1_fespace.GetVSize();
-      outsize = nd_fespace.GetVSize();
-   }
-   Vector xv(insize);
-   Vector assembled_y(outsize);
-   Vector pa_y(outsize);
-
-   xv.Randomize();
-   if (transpose)
-   {
-      assembled_grad_mat.BuildTranspose();
-      assembled_grad_mat.MultTranspose(xv, assembled_y);
-      pa_grad.MultTranspose(xv, pa_y);
-   }
-   else
-   {
-      assembled_grad_mat.Mult(xv, assembled_y);
-      pa_grad.Mult(xv, pa_y);
-   }
-
-   pa_y -= assembled_y;
-   double error = pa_y.Norml2() / assembled_y.Norml2();
-   INFO("dim " << dim << " ne " << num_elements << " order " << order
-        << (transpose ? " T:" : ":") << " error in PA gradient: " << error);
-
-   delete h1_fec;
-   delete nd_fec;
-   delete mesh;
-
-   return error;
-}
-
-TEST_CASE("PAGradient", "[CUDA]")
-{
-   auto transpose = GENERATE(true, false);
-   auto order = GENERATE(1, 2, 3, 4);
-   auto dim = GENERATE(2, 3);
-   auto num_elements = GENERATE(0, 1, 2, 3, 4);
-
-   double error = compare_pa_assembly(dim, num_elements, order, transpose);
-   REQUIRE(error == MFEM_Approx(0.0, 1.0e-14));
-}
-
-#ifdef MFEM_USE_MPI
-
-double par_compare_pa_assembly(int dim, int num_elements, int order,
-                               bool transpose)
-{
-   int rank;
-   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
-   int size;
-   MPI_Comm_size(MPI_COMM_WORLD, &size);
-
-   Mesh * smesh;
-   if (dim == 2)
-   {
-      smesh = new Mesh(num_elements, num_elements, Element::QUADRILATERAL, true);
-   }
-   else
-   {
-      smesh = new Mesh(num_elements, num_elements, num_elements,
-                       Element::HEXAHEDRON, true);
-   }
-   ParMesh * mesh = new ParMesh(MPI_COMM_WORLD, *smesh);
-   delete smesh;
-   FiniteElementCollection *h1_fec = new H1_FECollection(order, dim);
-   FiniteElementCollection *nd_fec = new ND_FECollection(order, dim);
-   ParFiniteElementSpace h1_fespace(mesh, h1_fec);
-   ParFiniteElementSpace nd_fespace(mesh, nd_fec);
-
-   ParDiscreteLinearOperator assembled_grad(&h1_fespace, &nd_fespace);
-   assembled_grad.AddDomainInterpolator(new GradientInterpolator);
-   const int skip_zeros = 1;
-   assembled_grad.Assemble(skip_zeros);
-   assembled_grad.Finalize(skip_zeros);
-   HypreParMatrix * assembled_grad_mat = assembled_grad.ParallelAssemble();
-
-   ParDiscreteLinearOperator pa_grad(&h1_fespace, &nd_fespace);
-   pa_grad.SetAssemblyLevel(AssemblyLevel::PARTIAL);
-   pa_grad.AddDomainInterpolator(new GradientInterpolator);
-   pa_grad.Assemble();
-   OperatorPtr pa_grad_oper;
-   pa_grad.FormRectangularSystemMatrix(pa_grad_oper);
-
-   int insize, outsize;
-   if (transpose)
-   {
-      insize = assembled_grad_mat->Height();
-      outsize = assembled_grad_mat->Width();
-   }
-   else
-   {
-      insize = assembled_grad_mat->Width();
-      outsize = assembled_grad_mat->Height();
-   }
-   Vector xv(insize);
-   Vector assembled_y(outsize);
-   Vector pa_y(outsize);
-   assembled_y = 0.0;
-   pa_y = 0.0;
-
-   xv.Randomize();
-   if (transpose)
-   {
-      assembled_grad_mat->MultTranspose(xv, assembled_y);
-      pa_grad_oper->MultTranspose(xv, pa_y);
-   }
-   else
-   {
-      assembled_grad_mat->Mult(xv, assembled_y);
-      pa_grad_oper->Mult(xv, pa_y);
-   }
-
-   Vector error_vec(pa_y);
-   error_vec -= assembled_y;
-   // serial norms and serial error; we are enforcing equality on each processor
-   // in the test
-   double error = error_vec.Norml2() / assembled_y.Norml2();
-
-   for (int p = 0; p < size; ++p)
-   {
-      if (rank == p)
-      {
-         INFO("[" << rank << "][par] dim " << dim << " ne " << num_elements
-              << " order " << order << (transpose ? " T:" : ":")
-              << " error in PA gradient: " << error);
-      }
-      MPI_Barrier(MPI_COMM_WORLD);
-   }
-
-   delete h1_fec;
-   delete nd_fec;
-   delete assembled_grad_mat;
-   delete mesh;
-
-   return error;
-}
-
-TEST_CASE("ParallelPAGradient", "[Parallel], [ParallelPAGradient]")
-{
-   auto transpose = GENERATE(true, false);
-   auto order = GENERATE(1, 2, 3, 4);
-   auto dim = GENERATE(2, 3);
-   auto num_elements = GENERATE(4, 5);
-
-   double error = par_compare_pa_assembly(dim, num_elements, order, transpose);
-   REQUIRE(error == MFEM_Approx(0.0, 1.0e-14));
-}
-
-#endif

tests/unit/fem/test_pa_idinterp.cpp

@@ -1,123 +0,0 @@
-// 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 "catch.hpp"
-#include "mfem.hpp"
-#include "unit_tests.hpp"
-
-using namespace mfem;
-
-double compare_pa_id_assembly(int dim, int num_elements, int order,
-                              bool transpose)
-{
-   Mesh * mesh;
-   if (num_elements == 0)
-   {
-      if (dim == 2)
-      {
-         mesh = new Mesh("../../data/star.mesh", order);
-      }
-      else
-      {
-         mesh = new Mesh("../../data/beam-hex.mesh", order);
-
-         // Transform mesh vertices to test without alignment with coordinate axes.
-         for (int i=0; i<mesh->GetNV(); ++i)
-         {
-            double *v = mesh->GetVertex(i);
-            const double yscale = 1.0 + v[1];
-            const double zscale = 1.0 + v[2];
-            v[0] *= zscale;
-            v[1] *= zscale;
-            v[2] *= yscale;
-         }
-      }
-   }
-   else
-   {
-      if (dim == 2)
-      {
-         mesh = new Mesh(num_elements, num_elements, Element::QUADRILATERAL, true);
-      }
-      else
-      {
-         mesh = new Mesh(num_elements, num_elements, num_elements,
-                         Element::HEXAHEDRON, true);
-      }
-   }
-   FiniteElementCollection *h1_fec = new H1_FECollection(order, dim);
-   FiniteElementCollection *nd_fec = new ND_FECollection(order, dim);
-   FiniteElementSpace h1_fespace(mesh, h1_fec, dim);
-   FiniteElementSpace nd_fespace(mesh, nd_fec);
-
-   DiscreteLinearOperator assembled_id(&h1_fespace, &nd_fespace);
-   assembled_id.AddDomainInterpolator(new IdentityInterpolator);
-   const int skip_zeros = 1;
-   assembled_id.Assemble(skip_zeros);
-   assembled_id.Finalize(skip_zeros);
-   const SparseMatrix& assembled_id_mat = assembled_id.SpMat();
-
-   DiscreteLinearOperator pa_id(&h1_fespace, &nd_fespace);
-   pa_id.SetAssemblyLevel(AssemblyLevel::PARTIAL);
-   pa_id.AddDomainInterpolator(new IdentityInterpolator);
-   pa_id.Assemble();
-   pa_id.Finalize();
-
-   int insize, outsize;
-   if (transpose)
-   {
-      insize = nd_fespace.GetVSize();
-      outsize = h1_fespace.GetVSize();
-   }
-   else
-   {
-      insize = h1_fespace.GetVSize();
-      outsize = nd_fespace.GetVSize();
-   }
-   Vector x(insize);
-   Vector assembled_y(outsize);
-   Vector pa_y(outsize);
-
-   x.Randomize();
-   if (transpose)
-   {
-      assembled_id_mat.BuildTranspose();
-      assembled_id_mat.MultTranspose(x, assembled_y);
-      pa_id.MultTranspose(x, pa_y);
-   }
-   else
-   {
-      assembled_id.Mult(x, assembled_y);
-      pa_id.Mult(x, pa_y);
-   }
-
-   pa_y -= assembled_y;
-   double error = pa_y.Norml2() / assembled_y.Norml2();
-   INFO("dim " << dim << " ne " << num_elements << " order " << order
-        << (transpose ? " T:" : ":") << " error in PA identity: " << error);
-
-   delete h1_fec;
-   delete nd_fec;
-   delete mesh;
-
-   return error;
-}
-
-TEST_CASE("PAIdentityInterp", "[CUDA]")
-{
-   auto transpose = GENERATE(true, false);
-   auto order = GENERATE(1, 2, 3, 4);
-   auto dim = GENERATE(2, 3);
-   auto num_elements = GENERATE(0, 1, 2, 3, 4);
-
-   double error = compare_pa_id_assembly(dim, num_elements, order, transpose);
-   REQUIRE(error == MFEM_Approx(0.0, 1.0e-14));
-}

tests/unit/fem/test_sum_bilin.cpp

@@ -1,231 +0,0 @@
-// 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