modifications

examples/CMakeLists.txt

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

examples/ex91.cpp

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

fem/CMakeLists.txt

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

fem/fem.hpp

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

fem/nonlininteg.cpp

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

fem/nonlininteg.hpp

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

fem/pprmnonlinearform.cpp

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

fem/pprmnonlinearform.hpp

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

fem/prmnonlinearform.hpp

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