oioi

miniapps/mtop/CMakeLists.txt

@@ -10,12 +10,10 @@
 # CONTRIBUTING.md for details.
 
 list(APPEND SEQMTOP_COMMON_SOURCES
-	paramnonlinearform.cpp
-        mtop_integrators.cpp)
+        )
 
 list(APPEND SEQMTOP_COMMON_HEADERS
-	paramnonlinearform.hpp
-        mtop_integrators.hpp)
+        )
 
 convert_filenames_to_full_paths(SEQMTOP_COMMON_SOURCES)
 convert_filenames_to_full_paths(SEQMTOP_COMMON_HEADERS)
@@ -24,17 +22,17 @@ set(SEQMTOP_COMMON_FILES
 	EXTRA_SOURCES ${SEQMTOP_COMMON_SOURCES}
 	EXTRA_HEADERS ${SEQMTOP_COMMON_HEADERS})
 
-add_mfem_miniapp(seqheat
-	MAIN seqheat.cpp
-	${SEQMTOP_COMMON_FILES}
-	LIBRARIES mfem)
 
 if(MFEM_USE_MPI)
 
 list(APPEND PARMTOP_COMMON_SOURCES
-	pparamnonlinearform.cpp)
+	rand_eigensolver.cpp
+        mtop_solvers.cpp
+        stokes_solver.cpp)
 list(APPEND PARMTOP_COMMON_HEADERS
-	pparamnonlinearform.hpp)
+	rand_eigensolver.hpp
+        mtop_solvers.hpp
+        stokes_solver.hpp)
 
 convert_filenames_to_full_paths(PARMTOP_COMMON_SOURCES)
 convert_filenames_to_full_paths(PARMTOP_COMMON_HEADERS)
@@ -46,8 +44,29 @@ set(PARMTOP_COMMON_FILES
 # message(STATUS "PARMTOP_COMMON_FILES: ${PARMTOP_COMMON_FILES}")
 # message(STATUS "SEQMTOP_COMMON_FILES: ${SEQMTOP_COMMON_FILES}")
 
-add_mfem_miniapp(parheat
-     MAIN parheat.cpp
+
+add_mfem_miniapp(test_eigp
+     MAIN test_eig.cpp
+     ${PARMTOP_COMMON_FILES}
+     LIBRARIES mfem)
+
+add_mfem_miniapp(test_stokes
+     MAIN test_stokes.cpp
+     ${PARMTOP_COMMON_FILES}
+     LIBRARIES mfem)
+
+add_mfem_miniapp(test_presb
+     MAIN test_presb.cpp
+     ${PARMTOP_COMMON_FILES}
+     LIBRARIES mfem)
+
+add_mfem_miniapp(test_msp1
+     MAIN test_msp1.cpp
+     ${PARMTOP_COMMON_FILES}
+     LIBRARIES mfem)
+
+add_mfem_miniapp(test_evec
+     MAIN test_evec.cpp
      ${PARMTOP_COMMON_FILES}
      LIBRARIES mfem)
 

miniapps/mtop/mtop_integrators.cpp

@@ -1,390 +0,0 @@
-// Copyright (c) 2010-2025, 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 "mtop_integrators.hpp"
-
-namespace mfem
-{
-
-real_t ParametricLinearDiffusion::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 dim = el[0]->GetDim();
-   int spaceDim = Tr.GetSpaceDim();
-   if (dim != spaceDim)
-   {
-      mfem::mfem_error("ParametricLinearDiffusion::GetElementEnergy"
-                       " is not defined on manifold meshes");
-   }
-
-   // shape functions
-   Vector shu0(dof_u0);
-   Vector shr0(dof_r0);
-   DenseMatrix dsu0(dof_u0,dim);
-   DenseMatrix B(dof_u0, 4);
-   B=0.0;
-
-   real_t w;
-
-   Vector param(1); param=0.0;
-   Vector uu(4); uu=0.0;
-
-   real_t energy =0.0;
-
-   const IntegrationRule *ir;
-   {
-      int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
-                 +pel[0]->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);
-
-      param[0]=shr0*(*pelfun[0]);
-
-      // 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 * qfun.QEnergy(Tr,ip,param,uu);
-   }
-   return energy;
-}
-
-
-void ParametricLinearDiffusion::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 dim = el[0]->GetDim();
-
-   elvec[0]->SetSize(dof_u0);
-   *elvec[0]=0.0;
-   int spaceDim = Tr.GetSpaceDim();
-   if (dim != spaceDim)
-   {
-      mfem::mfem_error("ParametricLinearDiffusion::AssembleElementVector"
-                       " is not defined on manifold meshes");
-   }
-
-   // shape functions
-   Vector shu0(dof_u0);
-   Vector shr0(dof_r0);
-   DenseMatrix dsu0(dof_u0,dim);
-   DenseMatrix B(dof_u0, 4);
-   B=0.0;
-
-   real_t w;
-
-   Vector param(1); param=0.0;
-   Vector uu(4); uu=0.0;
-   Vector rr(4);
-   Vector lvec; lvec.SetSize(dof_u0);
-
-   const IntegrationRule *ir = nullptr;
-   int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
-              +pel[0]->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);
-
-      param[0]=shr0*(*pelfun[0]);
-
-      // 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);
-      qfun.QResidual(Tr,ip,param, uu, rr);
-
-      B.Mult(rr,lvec);
-      elvec[0]->Add(w,lvec);
-   }
-}
-
-void ParametricLinearDiffusion::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 dim = el[0]->GetDim();
-
-   DenseMatrix* K=elmats(0,0);
-   K->SetSize(dof_u0,dof_u0);
-   (*K)=0.0;
-
-   int spaceDim = Tr.GetSpaceDim();
-   if (dim != spaceDim)
-   {
-      mfem::mfem_error("ParametricLinearDiffusion::AssembleElementGrad"
-                       " is not defined on manifold meshes");
-   }
-
-   // shape functions
-   Vector shu0(dof_u0);
-   Vector shr0(dof_r0);
-   DenseMatrix dsu0(dof_u0,dim);
-   DenseMatrix B(dof_u0, 4);
-   DenseMatrix A(dof_u0, 4);
-   B=0.0;
-   real_t w;
-
-   Vector param(1); param=0.0;
-   Vector uu(4); uu=0.0;
-   DenseMatrix hh(4,4);
-   Vector lvec; lvec.SetSize(dof_u0);
-
-   const IntegrationRule *ir = nullptr;
-   int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
-              +pel[0]->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);
-
-      param[0]=shr0*(*pelfun[0]);
-
-      // 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);
-      qfun.QGradResidual(Tr,ip,param,uu,hh);
-      Mult(B,hh,A);
-      AddMult_a_ABt(w,A,B,*K);
-   }
-}
-
-void ParametricLinearDiffusion::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 dim = el[0]->GetDim();
-   Vector& e0 = *(elvec[0]);
-
-   e0.SetSize(dof_r0);
-   e0=0.0;
-
-   int spaceDim = Tr.GetSpaceDim();
-   if (dim != spaceDim)
-   {
-      mfem::mfem_error("ParametricLinearDiffusion::AssemblePrmElementVector"
-                       " is not defined on manifold meshes");
-   }
-
-   // shape functions
-   Vector shu0(dof_u0);
-   Vector shr0(dof_r0);
-   DenseMatrix dsu0(dof_u0,dim);
-   DenseMatrix B(dof_u0, 4);
-   B=0.0;
-
-   real_t w;
-
-   Vector param(1); param=0.0;
-   Vector uu(4); uu=0.0;
-   Vector aa(4); aa=0.0;
-   Vector rr(1);
-   Vector lvec0; lvec0.SetSize(dof_r0);
-
-   const IntegrationRule *ir;
-   {
-      int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
-                 +pel[0]->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);
-
-      param[0]=shr0*(*pelfun[0]);
-
-      // 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);
-
-      qfun.AQResidual(Tr, ip, param, uu, aa, rr);
-
-      lvec0=shr0;
-      lvec0*=rr[0];
-
-      e0.Add(w,lvec0);
-   }
-}
-
-real_t DiffusionObjIntegrator::GetElementEnergy(const
-                                                Array<const FiniteElement *> &el,
-                                                ElementTransformation &Tr,
-                                                const Array<const Vector *> &elfun)
-{
-   int dof_u0 = el[0]->GetDof();
-   int dim = el[0]->GetDim();
-   int spaceDim = Tr.GetSpaceDim();
-   if (dim != spaceDim)
-   {
-      mfem::mfem_error("DiffusionObjIntegrator::GetElementEnergy"
-                       " is not defined on manifold meshes");
-   }
-
-   // shape functions
-   Vector shu0(dof_u0);
-
-   real_t w;
-   real_t val;
-
-   real_t energy = 0.0;
-
-   const IntegrationRule *ir;
-   {
-      int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0]);
-      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]->CalcPhysShape(Tr,shu0);
-
-      val=shu0*(*elfun[0]);
-      energy=energy + w * val * val;
-   }
-   return 0.5*energy;
-}
-
-void DiffusionObjIntegrator::AssembleElementVector(const
-                                                   Array<const FiniteElement *> &el,
-                                                   ElementTransformation &Tr,
-                                                   const Array<const Vector *> &elfun,
-                                                   const Array<Vector *> &elvec)
-{
-   int dof_u0 = el[0]->GetDof();
-   int dim = el[0]->GetDim();
-   int spaceDim = Tr.GetSpaceDim();
-
-   elvec[0]->SetSize(dof_u0);
-   *elvec[0]=0.0;
-
-   if (dim != spaceDim)
-   {
-      mfem::mfem_error("DiffusionObjIntegrator::GetElementEnergy"
-                       " is not defined on manifold meshes");
-   }
-
-   // shape functions
-   Vector shu0(dof_u0);
-
-   real_t w;
-   real_t val;
-
-   const IntegrationRule *ir;
-   {
-      int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0]);
-      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]->CalcPhysShape(Tr,shu0);
-
-      val=shu0*(*elfun[0]);
-
-      elvec[0]->Add(w*val,shu0);
-   }
-}
-
-
-} // end mfem namespace

miniapps/mtop/mtop_integrators.hpp

@@ -1,233 +0,0 @@
-// Copyright (c) 2010-2025, 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 MTOPINTEGRATORS_HPP
-#define MTOPINTEGRATORS_HPP
-
-#include "mfem.hpp"
-#include "paramnonlinearform.hpp"
-
-#include <map>
-
-namespace mfem
-{
-
-/// Base class for representing function at integration points.
-class BaseQFunction
-{
-public:
-   virtual ~BaseQFunction() {}
-
-   /// Returns a user defined string identifying the function.
-   virtual std::string GetType()=0;
-
-   // Returns the energy at an integration point.
-   virtual
-   real_t QEnergy(ElementTransformation &T, const IntegrationPoint &ip,
-                  mfem::Vector &dd, mfem::Vector &uu)
-   {
-      return 0.0;
-   }
-
-   // Returns the residual at an integration point.
-   virtual
-   void QResidual(ElementTransformation &T, const IntegrationPoint &ip,
-                  mfem::Vector &dd, mfem::Vector &uu, mfem::Vector &rr)=0;
-
-   /// Returns the gradient of the residual at a integration point.
-   virtual
-   void QGradResidual(ElementTransformation &T, const IntegrationPoint &ip,
-                      mfem::Vector &dd, mfem::Vector &uu, mfem::DenseMatrix &hh)=0;
-
-   /// Returns the gradient of the residual with respect to the design
-   /// parameters, multiplied by the adjoint.
-   virtual
-   void AQResidual(ElementTransformation &T, const IntegrationPoint &ip,
-                   mfem::Vector &dd, mfem::Vector &uu,
-                   mfem::Vector &aa, mfem::Vector &rr)=0;
-
-};
-
-/* QLinearDiffusion implements methods for computing the energy, the residual,
- * gradient of the residual and the product of the adjoint fields with the
- * derivative of the residual with respect to the parameters. All computations
- * are performed at a integration point. Therefore the vectors (vv,uu,aa,rr ..)
- * hold the fields' values and the fields' derivatives at the integration
- * point. For example for a single scalar parametric field representing the
- * density in topology optimization the vector dd will have size one and the
- * element will be the density at the integration point. The map between state
- * and parameter is not fixed and depends on the implementation of the QFunction
- * class. */
-class QLinearDiffusion:public BaseQFunction
-{
-public:
-   QLinearDiffusion(mfem::Coefficient& diffco, mfem::Coefficient& hsrco,
-                    real_t pp=1.0, real_t minrho=1e-7, real_t betac=4.0, real_t etac=0.5):
-      diff(diffco),load(hsrco), powerc(pp), rhomin(minrho), beta(betac), eta(etac)
-   {
-
-   }
-
-   std::string GetType() override
-   {
-      return "QLinearDiffusion";
-   }
-
-   real_t QEnergy(ElementTransformation &T, const IntegrationPoint &ip,
-                  Vector &dd, Vector &uu) override
-   {
-      // dd[0] - density
-      // uu[0] - grad_x
-      // uu[1] - grad_y
-      // uu[2] - grad_z
-      // uu[3] - temperature/scalar field
-
-      real_t di=diff.Eval(T,ip);
-      real_t ll=load.Eval(T,ip);
-      // Computes the physical density using projection.
-      real_t rz=0.5+0.5*std::tanh(beta*(dd[0]-eta)); //projection
-      // Computes the diffusion coefficient at the integration point.
-      real_t fd=di*(std::pow(rz,powerc)+rhomin);
-      // Computes the sum of the energy and the product of the temperature and
-      // the external input at the integration point.
-      real_t rez = 0.5*(uu[0]*uu[0]+uu[1]*uu[1]+uu[2]*uu[2])*fd-uu[3]*ll;
-      return rez;
-   }
-
-   /// Returns the derivative of QEnergy with respect to the state vector uu.
-   void QResidual(ElementTransformation &T, const IntegrationPoint &ip,
-                  Vector &dd, Vector &uu, Vector &rr) override
-   {
-      real_t di=diff.Eval(T,ip);
-      real_t ll=load.Eval(T,ip);
-      real_t rz=0.5+0.5*std::tanh(beta*(dd[0]-eta));
-      real_t fd=di*(std::pow(rz,powerc)+rhomin);
-
-      rr[0]=uu[0]*fd;
-      rr[1]=uu[1]*fd;
-      rr[2]=uu[2]*fd;
-      rr[3]=-ll;
-   }
-
-
-   // Returns the derivative, with respect to the density, of the product of
-   // the adjoint field with the residual at the integration point ip.
-   void AQResidual(ElementTransformation &T, const IntegrationPoint &ip,
-                   Vector &dd, Vector &uu, Vector &aa, Vector &rr) override
-   {
-      real_t di=diff.Eval(T,ip);
-      real_t tt=std::tanh(beta*(dd[0]-eta));
-      real_t rz=0.5+0.5*tt;
-      real_t fd=di*powerc*std::pow(rz,powerc-1.0)*0.5*(1.0-tt*tt)*beta;
-
-      rr[0] = -(aa[0]*uu[0]+aa[1]*uu[1]+aa[2]*uu[2])*fd;
-   }
-
-   // Returns the gradient of the residual with respect to the state vector at
-   // the integration point ip.
-   void QGradResidual(ElementTransformation &T, const IntegrationPoint &ip,
-                      Vector &dd, Vector &uu, DenseMatrix &hh) override
-   {
-      real_t di=diff.Eval(T,ip);
-      real_t tt=std::tanh(beta*(dd[0]-eta));
-      real_t rz=0.5+0.5*tt;
-      real_t fd=di*(std::pow(rz,powerc)+rhomin);
-      hh=0.0;
-
-      hh(0,0)=fd;
-      hh(1,1)=fd;
-      hh(2,2)=fd;
-      hh(3,3)=0.0;
-   }
-
-private:
-   mfem::Coefficient& diff; //diffusion coefficient
-   mfem::Coefficient& load; //load coefficient
-   real_t powerc; //penalization coefficient
-   real_t rhomin; //lower bound for the density
-   real_t beta;   //controls the sharpness of the projection
-   real_t eta;    //projection threshold for tanh
-};
-
-/// Provides implementation of an integrator for linear diffusion with
-/// parametrization provided by a density field. The setup is standard for
-/// topology optimization problems.
-class ParametricLinearDiffusion: public ParametricBNLFormIntegrator
-{
-public:
-   ParametricLinearDiffusion(BaseQFunction& qfunm): qfun(qfunm)
-   {
-
-   }
-
-   /// Computes the local energy.
-   real_t GetElementEnergy(const Array<const FiniteElement *> &el,
-                           const Array<const FiniteElement *> &pel,
-                           ElementTransformation &Tr,
-                           const Array<const Vector *> &elfun,
-                           const Array<const Vector *> &pelfun) override;
-
-   /// Computes the element's residual.
-   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) override;
-
-   /// Computes the stiffness/tangent matrix.
-   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) override;
-
-   /// Computes the product of the adjoint solution and the derivative of the
-   /// residual with respect to the parametric fields.
-   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) override;
-private:
-   BaseQFunction& qfun;
-};
-
-
-/// Computes an example of nonlinear objective
-/// $\int \rm{field}*\rm{field}*\rm{weight})\rm{d}\Omega_e$.
-class DiffusionObjIntegrator:public BlockNonlinearFormIntegrator
-{
-public:
-
-   DiffusionObjIntegrator()
-   {
-
-   }
-
-   /// Returns the objective contribution at element level.
-   real_t GetElementEnergy(const Array<const FiniteElement *> &el,
-                           ElementTransformation &Tr,
-                           const Array<const Vector *> &elfun) override;
-
-   /// Returns the gradient of the objective contribution at element level.
-   void AssembleElementVector(const Array<const FiniteElement *> &el,
-                              ElementTransformation &Tr,
-                              const Array<const Vector *> &elfun,
-                              const Array<Vector *> &elvec) override;
-};
-
-}
-
-#endif

miniapps/mtop/mtop_solvers.cpp

@@ -0,0 +1,395 @@
+#include "mfem.hpp"
+#include "mtop_solvers.hpp"
+#include "general/forall.hpp"
+
+
+LElasticOperator::LElasticOperator(mfem::ParMesh* mesh_, int vorder)
+{
+    pmesh=mesh_;
+    const int dim=pmesh->Dimension();
+    vfec=new mfem::H1_FECollection(vorder,dim);
+    vfes=new mfem::ParFiniteElementSpace(pmesh,vfec,dim, mfem::Ordering::byVDIM);
+
+    fdisp.SetSpace(vfes); fdisp=0.0;
+    adisp.SetSpace(vfes); adisp=0.0;
+
+    sol.SetSize(vfes->GetTrueVSize()); sol=0.0;
+    rhs.SetSize(vfes->GetTrueVSize()); rhs=0.0;
+    adj.SetSize(vfes->GetTrueVSize()); adj=0.0;
+
+    SetLinearSolver();
+
+    prec=nullptr;
+    ls=nullptr;
+
+    mfem::Operator::width=vfes->GetTrueVSize();
+    mfem::Operator::height=vfes->GetTrueVSize();
+
+    lvforce=nullptr;
+    volforce=nullptr;
+
+    E=nullptr;
+    nu=nullptr;
+
+    lambda=nullptr;
+    mu=nullptr;
+
+    bf=nullptr;
+    lf=nullptr;
+}
+
+LElasticOperator::~LElasticOperator()
+{
+    delete prec;
+    delete ls;
+
+    delete bf;
+    delete lf;
+
+    delete vfes;
+    delete vfec;
+
+    delete lvforce;
+
+    for(auto it=load_coeff.begin();it!=load_coeff.end();it++){
+        delete it->second;
+    }
+
+    delete lambda;
+    delete mu;
+}
+
+void LElasticOperator::SetLinearSolver(mfem::real_t rtol, mfem::real_t atol, int miter)
+{
+    linear_rtol=rtol;
+    linear_atol=atol;
+    linear_iter=miter;
+}
+
+void LElasticOperator::AddDispBC(int id, int dir, mfem::real_t val)
+{
+    if(dir==0){
+        bcx[id]=mfem::ConstantCoefficient(val);
+        AddDispBC(id,dir,bcx[id]);
+    }
+    if(dir==1){
+        bcy[id]=mfem::ConstantCoefficient(val);
+        AddDispBC(id,dir,bcy[id]);
+
+    }
+    if(dir==2){
+        bcz[id]=mfem::ConstantCoefficient(val);
+        AddDispBC(id,dir,bcz[id]);
+    }
+    if(dir==4){
+        bcx[id]=mfem::ConstantCoefficient(val);
+        bcy[id]=mfem::ConstantCoefficient(val);
+        bcz[id]=mfem::ConstantCoefficient(val);
+        AddDispBC(id,0,bcx[id]);
+        AddDispBC(id,1,bcy[id]);
+        AddDispBC(id,2,bcz[id]);
+    }
+}
+
+void LElasticOperator::DelDispBC()
+{
+    bccx.clear();
+    bccy.clear();
+    bccz.clear();
+
+    bcx.clear();
+    bcy.clear();
+    bcz.clear();
+
+    ess_tdofv.DeleteAll();
+}
+
+void LElasticOperator::AddDispBC(int id, int dir, mfem::Coefficient &val)
+{
+    if(dir==0){ bccx[id]=&val; }
+    if(dir==1){ bccy[id]=&val; }
+    if(dir==2){ bccz[id]=&val; }
+    if(dir==4){ bccx[id]=&val; bccy[id]=&val; bccz[id]=&val;}
+    if(pmesh->Dimension()==2)
+    {
+        bccz.clear();
+    }
+}
+
+void LElasticOperator::SetVolForce(mfem::real_t fx, double fy, double fz)
+{
+    delete lvforce;
+    int dim=pmesh->Dimension();
+    mfem::Vector ff(dim); ff(0)=fx; ff(1)=fy;
+    if(dim==3){ff(2)=fz;}
+    lvforce=new mfem::VectorConstantCoefficient(ff);
+    volforce=lvforce;
+
+}
+
+void LElasticOperator::SetVolForce(mfem::VectorCoefficient& fv)
+{
+    volforce=&fv;
+}
+
+void LElasticOperator::SetEssTDofs(mfem::Vector& bsol, mfem::Array<int>& ess_dofs)
+{
+    // Set the BC
+
+    ess_tdofv.DeleteAll();
+
+    mfem::Array<int> ess_tdofx;
+    mfem::Array<int> ess_tdofy;
+    mfem::Array<int> ess_tdofz;
+
+    int dim=pmesh->Dimension();
+
+    {
+        for(auto it=bccx.begin();it!=bccx.end();it++)
+        {
+            mfem::Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+            ess_bdr=0;
+            ess_bdr[it->first -1]=1;
+            mfem::Array<int> ess_tdof_list;
+            vfes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list,0);
+            ess_tdofx.Append(ess_tdof_list);
+
+            mfem::VectorArrayCoefficient pcoeff(dim);
+            pcoeff.Set(0, it->second, false);
+            fdisp.ProjectBdrCoefficient(pcoeff, ess_bdr);
+        }
+        //copy tdofsx from velocity grid function
+        {
+            mfem::Vector& vc=fdisp.GetTrueVector();
+            for(int ii=0;ii<ess_tdofx.Size();ii++)
+            {
+                bsol[ess_tdofx[ii]]=vc[ess_tdofx[ii]];
+            }
+        }
+        ess_dofs.Append(ess_tdofx); ess_tdofx.DeleteAll();
+
+        for(auto it=bccy.begin();it!=bccy.end();it++)
+        {
+            mfem::Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+            ess_bdr=0;
+            ess_bdr[it->first -1]=1;
+            mfem::Array<int> ess_tdof_list;
+            vfes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list,1);
+            ess_tdofy.Append(ess_tdof_list);
+
+            mfem::VectorArrayCoefficient pcoeff(dim);
+            pcoeff.Set(1, it->second, false);
+            fdisp.ProjectBdrCoefficient(pcoeff, ess_bdr);
+        }
+        //copy tdofsy from velocity grid function
+        {
+            mfem::Vector& vc=fdisp.GetTrueVector();
+            for(int ii=0;ii<ess_tdofy.Size();ii++)
+            {
+                bsol[ess_tdofy[ii]]=vc[ess_tdofy[ii]];
+            }
+        }
+        ess_dofs.Append(ess_tdofy); ess_tdofy.DeleteAll();
+
+        if(dim==3){
+            for(auto it=bccz.begin();it!=bccz.end();it++)
+            {
+                mfem::Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+                ess_bdr=0;
+                ess_bdr[it->first -1]=1;
+                mfem::Array<int> ess_tdof_list;
+                vfes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list,2);
+                ess_tdofz.Append(ess_tdof_list);
+
+                mfem::VectorArrayCoefficient pcoeff(dim);
+                pcoeff.Set(2, it->second, false);
+                fdisp.ProjectBdrCoefficient(pcoeff, ess_bdr);
+            }
+
+            //copy tdofsz from velocity grid function
+            {
+                mfem::Vector& vc=fdisp.GetTrueVector();
+                for(int ii=0;ii<ess_tdofz.Size();ii++)
+                {
+                    bsol[ess_tdofz[ii]]=vc[ess_tdofz[ii]];
+                }
+            }
+            ess_dofs.Append(ess_tdofz); ess_tdofz.DeleteAll();
+        }
+    }
+
+}
+
+void LElasticOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const
+{
+    //the rhs x is assumed to have the contribution of the BC
+    ls->Mult(x,y);
+
+    int N=ess_tdofv.Size();
+    mfem::real_t *yp = y.ReadWrite();
+    const mfem::real_t *sp = sol.Read();
+    const int *ep = ess_tdofv.Read();
+    mfem::forall(N, [=] MFEM_HOST_DEVICE (int i)
+    {
+        yp[ep[i]]=sp[ep[i]];
+    });
+}
+
+void LElasticOperator::MultTranspose(const mfem::Vector &x, mfem::Vector &y) const
+{
+    //the adjoint rhs is assumed to be corrected for the BC
+    //K is symmetric
+    ls->Mult(x,y);
+
+    int N=ess_tdofv.Size();
+    ess_tdofv.Read();
+
+    mfem::real_t *yp = y.ReadWrite();
+    const int *ep = ess_tdofv.Read();
+
+    mfem::forall(N,[=] MFEM_HOST_DEVICE (int i)
+    {
+        yp[ep[i]]=mfem::real_t(0.0);
+    });
+}
+
+void LElasticOperator::Assemble()
+{
+    if(bf==nullptr){return;}
+
+    //set BC
+    sol=mfem::real_t(0.0);
+    SetEssTDofs(sol,ess_tdofv);
+
+    bf->Assemble(0);
+    bf->FormSystemMatrix(ess_tdofv,K);
+
+    if(ls==nullptr){
+        ls=new CGSolver(pmesh->GetComm());
+        ls->SetAbsTol(linear_atol);
+        ls->SetRelTol(linear_rtol);
+        ls->SetMaxIter(linear_iter);
+        prec=new mfem::HypreBoomerAMG();
+        //set the rigid body modes
+        prec->SetElasticityOptions(vfes);
+        prec->SetPrintLevel(1);
+        ls->SetPreconditioner(*prec);
+        ls->SetOperator(K);
+        ls->SetPrintLevel(1);
+    }else{
+        ls->SetOperator(K);
+    }
+
+    std::cout<<pmesh->GetMyRank()<<" LSW="<<ls->Width()<<" LSH="<<ls->Height()
+                                 <<" KFW="<<K.Width()<<" KFH="<<K.Height()<<std::endl;
+}
+
+void LElasticOperator::FSolve()
+{
+    if(lf==nullptr){
+        lf=new ParLinearForm(vfes);
+        if(volforce!=nullptr){
+            lf->AddDomainIntegrator(new VectorDomainLFIntegrator(*volforce));
+        }
+        //add surface loads
+
+    }
+
+    (*lf)=mfem::real_t(0.0);
+
+    lf->Assemble();
+    lf->ParallelAssemble(rhs);
+
+    for(int i=0;i<ess_tdofv.Size();i++){
+        rhs[ess_tdofv[i]]=sol[ess_tdofv[i]];
+    }
+
+    ls->Mult(rhs,sol);
+
+}
+
+FRElasticSolver::FRElasticSolver(mfem::ParMesh* mesh_, int vorder, mfem::real_t freq_):
+                    LElasticOperator(mesh_,vorder)
+{
+
+    freq=freq_;
+    alpha=0.0;
+    beta=0.0;
+
+    mf=nullptr;
+    cf=nullptr;
+
+    lrf=nullptr;
+    lif=nullptr;
+
+    num_svd_modes=5;
+    num_svd_iter=6;
+    pop=nullptr;
+    ss_solver=nullptr;
+}
+
+FRElasticSolver::~FRElasticSolver()
+{
+
+    delete ss_solver;
+    delete pop;
+
+    delete mf;
+    delete cf;
+
+    delete lrf;
+    delete lif;
+
+}
+
+void FRElasticSolver::Assemble()
+{
+    LElasticOperator::Assemble();
+
+    if(mf==nullptr){
+        MFEM_WARNING("FRElasticSolver::Mass bilinear form is not defined!");
+        return;
+    }
+    mf->Assemble();
+    mf->FormSystemMatrix(LElasticOperator::ess_tdofv,hmf);
+
+    if(cf==nullptr){
+        MFEM_WARNING("FRElasticSolver::Damping bilinear form is not defined!");
+        return;
+    }
+    cf->Assemble();
+    cf->FormSystemMatrix(LElasticOperator::ess_tdofv,hcf);
+}
+
+void FRElasticSolver::AssembleSVD()
+{
+    delete ss_solver;
+    delete pop;
+    ss_solver=new RandomizedSubspaceIteration(LElasticOperator::pmesh->GetComm());
+
+    std::cout<<LElasticOperator::pmesh->GetMyRank()<<"CW="<<hmf->Width()<<" CH="<<hmf->Height()<<std::endl;
+
+    pop=new LocProductOperator(LElasticOperator::ls,hmf.Ptr());
+    ss_solver->SetOperator(*pop);
+    ss_solver->SetNumModes(num_svd_modes);
+    ss_solver->SetNumIter(num_svd_iter);
+    ss_solver->SetConstrDOFs(LElasticOperator::ess_tdofv);
+    ss_solver->Solve();
+
+    auto eigv=ss_solver->GetModes();
+    Vector tmpv; tmpv.SetSize(eigv[0].Size());
+    for(int i=0;i<ss_solver->GetNumModes();i++){
+        pop->MultB(eigv[i],tmpv);
+        for(int j=0;j<ss_solver->GetNumModes();j++){
+             real_t bb=InnerProduct(LElasticOperator::pmesh->GetComm(),tmpv,eigv[j]);
+             if(0==pmesh->GetMyRank()){
+                 std::cout<<bb<<" ";
+             }
+
+        }
+        if(0==pmesh->GetMyRank()){
+            std::cout<<std::endl;
+        }
+    }
+
+}

miniapps/mtop/paramnonlinearform.hpp

@@ -1,300 +0,0 @@
-// Copyright (c) 2010-2025, 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 "mfem.hpp"
-
-namespace mfem
-{
-
-/** The abstract base class ParametricBNLFormIntegrator is a generalization of
-    the BlockNonlinearFormIntegrator class suitable for block state and
-    parameter vectors. */
-class ParametricBNLFormIntegrator
-{
-public:
-   /// Compute the local energy
-   virtual real_t 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);
-
-   /// Perform the local action of the BlockNonlinearFormIntegrator on element
-   /// faces
-   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 BNLFormIntegrator
-   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 *> &pelvec);
-
-   /// Perform the local action on the parameters of the BNLFormIntegrator on
-   /// faces
-   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 *> &pelvect);
-
-   /// 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);
-
-   /// Assemble the local gradient matrix on faces of the elements
-   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 ~ParametricBNLFormIntegrator() { }
-};
-
-
-/** @brief A class representing a general parametric block nonlinear operator
-    defined on the Cartesian product of multiple FiniteElementSpace%s. */
-class ParametricBNLForm : public Operator
-{
-protected:
-   /// FE spaces on which the form lives.
-   Array<FiniteElementSpace*> fes;
-
-   /// FE spaces for the parametric fields
-   Array<FiniteElementSpace*> paramfes;
-
-   int paramheight;
-   int paramwidth;
-
-   /// Set of Domain Integrators to be assembled (added).
-   Array<ParametricBNLFormIntegrator*> dnfi;
-
-   /// Set of interior face Integrators to be assembled (added).
-   Array<ParametricBNLFormIntegrator*> fnfi;
-
-   /// Set of Boundary Face Integrators to be assembled (added).
-   Array<ParametricBNLFormIntegrator*> 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> paramblock_offsets;
-   Array<int> paramblock_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 'paramfes'
-   Array<Array<int> *> paramess_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 paramfes, may be NULL
-   Array<const Operator *> Pparam;
-
-   /// Array of results of dynamic-casting P to SparseMatrix pointer
-   Array<const SparseMatrix *> cP;
-
-   /// Array of results of dynamic-casting Pparam to SparseMatrix pointer
-   Array<const SparseMatrix *> cPparam;
-
-   /// 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 &ParamProlongate(const BlockVector &bx) const;
-
-   real_t 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 MultParamBlocked(const BlockVector &bx, const BlockVector & ax,
-                         const BlockVector &dx, BlockVector &dy) const;
-
-
-   /// Specialized version of GetGradient() for BlockVector
-   void ComputeGradientBlocked(const BlockVector &bx, const BlockVector &dx) const;
-
-public:
-   /// Construct an empty BlockNonlinearForm. Initialize with SetSpaces().
-   ParametricBNLForm();
-
-   /// Construct a BlockNonlinearForm on the given set of FiniteElementSpace%s.
-   ParametricBNLForm(Array<FiniteElementSpace *> &statef,
-                     Array<FiniteElementSpace *> &paramf);
-
-   /// Return the @a k-th FE space of the ParametricBNLForm.
-   FiniteElementSpace *FESpace(int k) { return fes[k]; }
-
-   /// Return the @a k-th parametric FE space of the ParametricBNLForm.
-   FiniteElementSpace *ParamFESpace(int k) { return paramfes[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 *ParamFESpace(int k) const { return paramfes[k]; }
-
-   /// Return the integrators
-   Array<ParametricBNLFormIntegrator*>& GetDNFI() { return dnfi;}
-
-
-   /// (Re)initialize the ParametricBNLForm.
-   /** After a call to SetSpaces(), the essential b.c. must be set again. */
-   void SetSpaces(Array<FiniteElementSpace *> &statef,
-                  Array<FiniteElementSpace *> &paramf);
-
-   /// 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> &ParamGetBlockOffsets() const { return paramblock_offsets; }
-
-   /// Return the true-dof offsets for the parameters.
-   const Array<int> &ParamGetBlockTrueOffsets() const { return paramblock_trueOffsets; }
-
-   /// Adds new Domain Integrator.
-   void AddDomainIntegrator(ParametricBNLFormIntegrator *nlfi)
-   { dnfi.Append(nlfi); }
-
-   /// Adds new Interior Face Integrator.
-   void AddInteriorFaceIntegrator(ParametricBNLFormIntegrator *nlfi)
-   { fnfi.Append(nlfi); }
-
-   /// Adds new Boundary Face Integrator.
-   void AddBdrFaceIntegrator(ParametricBNLFormIntegrator *nlfi)
-   { bfnfi.Append(nlfi); bfnfi_marker.Append(NULL); }
-
-   /** @brief Adds new Boundary Face Integrator, restricted to specific boundary
-       attributes. */
-   void AddBdrFaceIntegrator(ParametricBNLFormIntegrator *nlfi,
-                             Array<int> &bdr_marker);
-
-   /// Set the essential boundary conditions.
-   virtual void SetEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
-                               Array<Vector *> &rhs);
-
-   /// Set the essential boundary conditions on the parametric fields.
-   virtual void SetParamEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
-                                    Array<Vector *> &rhs);
-
-
-   /// Computes the energy for a state vector x.
-   virtual real_t GetEnergy(const Vector &x) const;
-
-   /// Method is only called in serial, the parallel version calls MultBlocked
-   /// directly.
-   void Mult(const Vector &x, Vector &y) const override;
-
-   /// Method is only called in serial, the parallel version calls MultBlocked
-   /// directly.
-   virtual void ParamMult(const Vector &x, Vector &y) const;
-
-   /// Method is only called in serial, the parallel version calls
-   /// GetGradientBlocked directly.
-   BlockOperator &GetGradient(const Vector &x) const override;
-
-   /// 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 SetParamFields(const Vector &dv) const;
-
-   /// Destructor.
-   virtual ~ParametricBNLForm();
-
-};
-
-}
-
-#endif

miniapps/mtop/parheat.cpp

@@ -1,354 +0,0 @@
-// Copyright (c) 2010-2025, 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.
-//
-//     ----------------------------------------------------------------
-//     ParHeat Miniapp: Gradients of PDE constrained objective function
-//     ----------------------------------------------------------------
-//                           (Parallel Version)
-//
-// The following example computes the gradients of a specified objective
-// function with respect to parametric fields. The objective function is having
-// the following form f(u(\rho)) where u(\rho) is a solution of a specific state
-// problem (in the example that is the diffusion equation), and \rho is a
-// parametric field discretized by finite elements. The parametric field (also
-// called density in topology optimization) controls the coefficients of the
-// state equation. For the considered case, the density controls the diffusion
-// coefficient within the computational domain.
-//
-// For more information, the users are referred to:
-//
-//    Hinze, M.; Pinnau, R.; Ulbrich, M. & Ulbrich, S.
-//    Optimization with PDE Constraints
-//    Springer Netherlands, 2009
-//
-//    Bendsøe, M. P. & Sigmund, O.
-//    Topology Optimization - Theory, Methods and Applications
-//    Springer Verlag, Berlin Heidelberg, 2003
-//
-// Compile with: make parheat
-//
-// Sample runs:
-//
-//    mpirun -np 4 parheat --visualization
-//    mpirun -np 4 parheat --visualization -m ../../data/beam-quad.mesh
-
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-
-#include "pparamnonlinearform.hpp"
-#include "mtop_integrators.hpp"
-
-using namespace mfem;
-
-int main(int argc, char *argv[])
-{
-   // 1. Initialize MPI and HYPRE.
-   mfem::Mpi::Init(argc, argv);
-   int myrank = mfem::Mpi::WorldRank();
-   mfem::Hypre::Init();
-
-   // Parse command-line options.
-   const char *mesh_file = "../../data/star.mesh";
-   int order = 1;
-   bool static_cond = false;
-   int ser_ref_levels = 1;
-   int par_ref_levels = 1;
-   real_t newton_rel_tol = 1e-7;
-   real_t newton_abs_tol = 1e-12;
-   int newton_iter = 10;
-   int print_level = 1;
-   bool visualization = false;
-
-   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(&par_ref_levels,
-                  "-rp",
-                  "--refine-parallel",
-                  "Number of times to refine the mesh uniformly in parallel.");
-   args.AddOption(&order, "-o", "--order",
-                  "Finite element order (polynomial degree) or -1 for"
-                  " isoparametric space.");
-   args.AddOption(&visualization,
-                  "-vis",
-                  "--visualization",
-                  "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
-                  "--no-static-condensation", "Enable static condensation.");
-   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())
-   {
-      if (myrank == 0)
-      {
-         args.PrintUsage(std::cout);
-      }
-      return 1;
-   }
-
-   if (myrank == 0)
-   {
-      args.PrintOptions(std::cout);
-   }
-
-   // Read the (serial) mesh from the given mesh file on all processors.  We
-   // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
-   // and volume meshes with the same code.
-   mfem::Mesh mesh(mesh_file, 1, 1);
-   int dim = mesh.Dimension();
-
-   // Refine the serial mesh on all processors to increase the resolution. In
-   // this example we do 'ref_levels' of uniform refinement. We choose
-   // 'ref_levels' to be the largest number that gives a final mesh with no
-   // more than 10,000 elements.
-   {
-      int ref_levels =
-         (int)floor(log(10000./mesh.GetNE())/log(2.)/dim);
-      for (int l = 0; l < ref_levels; l++)
-      {
-         mesh.UniformRefinement();
-      }
-   }
-
-   // Define a parallel mesh by a partitioning of the serial mesh. Refine
-   // this mesh further in parallel to increase the resolution. Once the
-   // parallel mesh is defined, the serial mesh can be deleted.
-   mfem::ParMesh pmesh(MPI_COMM_WORLD, mesh);
-   mesh.Clear();
-   {
-      for (int l = 0; l < par_ref_levels; l++)
-      {
-         pmesh.UniformRefinement();
-      }
-   }
-
-   // Define the Diffusion coefficient.
-   mfem::ConstantCoefficient* diffco=new mfem::ConstantCoefficient(1.0);
-   // Define the Heat source.
-   mfem::ConstantCoefficient* loadco=new mfem::ConstantCoefficient(1.0);
-   // Define the q-function.
-   mfem::QLinearDiffusion* qfun=new mfem::QLinearDiffusion(*diffco,*loadco,1.0,
-                                                           1e-7,4.0,0.5);
-
-   // Define FE collection and space for the state solution.
-   mfem::H1_FECollection sfec(order, dim);
-   mfem::ParFiniteElementSpace* sfes=new mfem::ParFiniteElementSpace(&pmesh,&sfec,
-                                                                     1);
-   // Define FE collection and space for the density field.
-   mfem::L2_FECollection pfec(order, dim);
-   mfem::ParFiniteElementSpace* pfes=new mfem::ParFiniteElementSpace(&pmesh,&pfec,
-                                                                     1);
-
-   // Define the arrays for the nonlinear form.
-   mfem::Array<mfem::ParFiniteElementSpace*> asfes;
-   mfem::Array<mfem::ParFiniteElementSpace*> apfes;
-
-   asfes.Append(sfes);
-   apfes.Append(pfes);
-
-   // Define parametric block nonlinear form using single scalar H1 field
-   // and L2 scalar density field.
-   mfem::ParParametricBNLForm* nf=new mfem::ParParametricBNLForm(asfes,apfes);
-   // Add a parametric integrator.
-   nf->AddDomainIntegrator(new mfem::ParametricLinearDiffusion(*qfun));
-
-   // Define true block vectors for state, adjoint, resudual.
-   mfem::BlockVector solbv; solbv.Update(nf->GetBlockTrueOffsets());    solbv=0.0;
-   mfem::BlockVector adjbv; adjbv.Update(nf->GetBlockTrueOffsets());    adjbv=0.0;
-   mfem::BlockVector resbv; resbv.Update(nf->GetBlockTrueOffsets());    resbv=0.0;
-   // Define true block vectors for parametric field and gradients.
-   mfem::BlockVector prmbv; prmbv.Update(nf->ParamGetBlockTrueOffsets());
-   prmbv=0.0;
-   mfem::BlockVector grdbv; grdbv.Update(nf->ParamGetBlockTrueOffsets());
-   grdbv=0.0;
-
-   // Set the BCs for the physics.
-   mfem::Array<mfem::Array<int> *> ess_bdr;
-   mfem::Array<mfem::Vector*>      ess_rhs;
-   ess_bdr.Append(new mfem::Array<int>(pmesh.bdr_attributes.Max()));
-   ess_rhs.Append(nullptr);
-   (*ess_bdr[0]) = 1;
-   nf->SetEssentialBC(ess_bdr,ess_rhs);
-   delete ess_bdr[0];
-
-   // Set the density field to 0.5.
-   prmbv=0.5;
-   // Set the density as parametric field in the parametric BNLForm.
-   nf->SetParamFields(prmbv); //set the density
-
-   // Compute the stiffness/tangent matrix for density prmbv=0.5.
-   mfem::BlockOperator *A = &nf->GetGradient(solbv);
-   mfem::HypreBoomerAMG* prec=new mfem::HypreBoomerAMG();
-   prec->SetPrintLevel(print_level);
-   // Use only block (0,0) as in this case we have a single field.
-   prec->SetOperator(A->GetBlock(0,0));
-
-   // Construct block preconditioner for the BNLForm.
-   mfem::BlockDiagonalPreconditioner *blpr = new mfem::BlockDiagonalPreconditioner(
-      nf->GetBlockTrueOffsets());
-   blpr->SetDiagonalBlock(0,prec);
-
-   // Define the solvers.
-   mfem::GMRESSolver *gmres;
-   gmres = new mfem::GMRESSolver(MPI_COMM_WORLD);
-   gmres->SetAbsTol(newton_abs_tol/10);
-   gmres->SetRelTol(newton_rel_tol/10);
-   gmres->SetMaxIter(100);
-   gmres->SetPrintLevel(print_level);
-   gmres->SetPreconditioner(*blpr);
-   gmres->SetOperator(*A);
-
-
-   // Solve the problem.
-   solbv=0.0;
-   nf->Mult(solbv,resbv); resbv.Neg(); //compute RHS
-   gmres->Mult(resbv, solbv);
-
-   // Compute the energy of the state system.
-   real_t energy = nf->GetEnergy(solbv);
-   if (myrank==0)
-   {
-      std::cout << "energy =" << energy << std::endl;
-   }
-
-   // Define the block nonlinear form utilized for representing the objective -
-   // use the state array from the BNLForm.
-   mfem::ParBlockNonlinearForm* ob=new mfem::ParBlockNonlinearForm(asfes);
-   // Add the integrator for the objective.
-   ob->AddDomainIntegrator(new mfem::DiffusionObjIntegrator());
-
-   // Compute the objective.
-   real_t obj=ob->GetEnergy(solbv);
-   if (myrank==0)
-   {
-      std::cout << "Objective =" << obj << std::endl;
-   }
-
-   // Solve the adjoint.
-   {
-      mfem::BlockVector adjrhs; adjrhs.Update(nf->GetBlockTrueOffsets());  adjrhs=0.0;
-      // Compute the RHS for the adjoint, i.e., the gradients with respect to
-      // the parametric fields.
-      ob->Mult(solbv, adjrhs);
-      // Get the tangent matrix from the state problem. We do not need to
-      // transpose the operator for diffusion.  Compute the adjoint solution.
-      gmres->Mult(adjrhs, adjbv);
-   }
-
-   // Compute gradients.
-   // First set the adjoint field.
-   nf->SetAdjointFields(adjbv);
-   // Set the state field.
-   nf->SetStateFields(solbv);
-   // Call the parametric Mult.
-   nf->ParamMult(prmbv, grdbv);
-
-   // Dump out the data.
-   if (visualization)
-   {
-      mfem::ParaViewDataCollection *dacol=new mfem::ParaViewDataCollection("ParHeat",
-                                                                           &pmesh);
-      mfem::ParGridFunction gfgrd(pfes); gfgrd.SetFromTrueDofs(grdbv.GetBlock(0));
-      mfem::ParGridFunction gfdns(pfes); gfdns.SetFromTrueDofs(prmbv.GetBlock(0));
-      // Define state grid function.
-      mfem::ParGridFunction gfsol(sfes); gfsol.SetFromTrueDofs(solbv.GetBlock(0));
-      mfem::ParGridFunction gfadj(sfes); gfadj.SetFromTrueDofs(adjbv.GetBlock(0));
-
-      dacol->SetLevelsOfDetail(order);
-      dacol->RegisterField("sol", &gfsol);
-      dacol->RegisterField("adj", &gfadj);
-      dacol->RegisterField("dns", &gfdns);
-      dacol->RegisterField("grd", &gfgrd);
-
-      dacol->SetTime(1.0);
-      dacol->SetCycle(1);
-      dacol->Save();
-
-      delete dacol;
-   }
-
-   // FD check
-   {
-      mfem::BlockVector prtbv;
-      mfem::BlockVector tmpbv;
-      prtbv.Update(nf->ParamGetBlockTrueOffsets());
-      tmpbv.Update(nf->ParamGetBlockTrueOffsets());
-      prtbv.GetBlock(0).Randomize();
-      prtbv*=1.0;
-      real_t lsc=1.0;
-
-      real_t gQoI=ob->GetEnergy(solbv);
-      real_t lQoI;
-
-      real_t nd=mfem::InnerProduct(MPI_COMM_WORLD,prtbv,prtbv);
-      real_t td=mfem::InnerProduct(MPI_COMM_WORLD,prtbv,grdbv);
-      td=td/nd;
-
-      for (int l = 0; l < 10; l++)
-      {
-         lsc/=10.0;
-         prtbv/=10.0;
-         add(prmbv,prtbv,tmpbv);
-         nf->SetParamFields(tmpbv);
-         // Solve the physics.
-         solbv=0.0;
-         nf->Mult(solbv,resbv); resbv.Neg(); //compute RHS
-         A = &nf->GetGradient(solbv);
-         prec->SetPrintLevel(0);
-         prec->SetOperator(A->GetBlock(0,0));
-         gmres->SetOperator(*A);
-         gmres->SetPrintLevel(0);
-         gmres->Mult(resbv,solbv);
-         // Compute the objective.
-         lQoI=ob->GetEnergy(solbv);
-         real_t ld=(lQoI-gQoI)/lsc;
-         if (myrank==0)
-         {
-            std::cout << "dx=" << lsc <<" FD approximation=" << ld/nd
-                      << " adjoint gradient=" << td
-                      << " err=" << std::fabs(ld/nd-td) << std::endl;
-         }
-      }
-   }
-
-   delete ob;
-   delete gmres;
-   delete blpr;
-   delete prec;
-
-   delete nf;
-   delete pfes;
-   delete sfes;
-
-   delete qfun;
-   delete loadco;
-   delete diffco;
-
-   return 0;
-}

miniapps/mtop/pparamnonlinearform.cpp

@@ -1,362 +0,0 @@
-// Copyright (c) 2010-2025, 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 "pparamnonlinearform.hpp"
-
-#ifdef MFEM_USE_MPI
-
-namespace mfem
-{
-
-ParParametricBNLForm::ParParametricBNLForm(Array<ParFiniteElementSpace *>
-                                           &statef,
-                                           Array<ParFiniteElementSpace *> &paramf)
-   :ParametricBNLForm()
-{
-   pBlockGrad = nullptr;
-   SetParSpaces(statef,paramf);
-}
-
-void ParParametricBNLForm::SetParSpaces(Array<ParFiniteElementSpace *> &statef,
-                                        Array<ParFiniteElementSpace *> &paramf)
-{
-   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(statef.Size());
-   Array<FiniteElementSpace *> prmserialSpaces(paramf.Size());
-   for (int s=0; s<statef.Size(); s++)
-   {
-      serialSpaces[s] = (FiniteElementSpace *) statef[s];
-   }
-   for (int s=0; s<paramf.Size(); s++)
-   {
-      prmserialSpaces[s] = (FiniteElementSpace *) paramf[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 * ParParametricBNLForm::ParFESpace(int k)
-{
-   return (ParFiniteElementSpace *)fes[k];
-}
-
-const ParFiniteElementSpace *ParParametricBNLForm::ParFESpace(int k) const
-{
-   return (const ParFiniteElementSpace *)fes[k];
-}
-
-
-ParFiniteElementSpace * ParParametricBNLForm::ParParamFESpace(int k)
-{
-   return (ParFiniteElementSpace *)paramfes[k];
-}
-
-const ParFiniteElementSpace *ParParametricBNLForm::ParParamFESpace(int k) const
-{
-   return (const ParFiniteElementSpace *)paramfes[k];
-}
-
-// Here, rhs is a true dof vector
-void ParParametricBNLForm::SetEssentialBC(const
-                                          Array<Array<int> *>&bdr_attr_is_ess,
-                                          Array<Vector *> &rhs)
-{
-   Array<Vector *> nullarray(fes.Size());
-   nullarray = NULL;
-
-   ParametricBNLForm::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 ParParametricBNLForm::SetParamEssentialBC(const
-                                               Array<Array<int> *>&bdr_attr_is_ess,
-                                               Array<Vector *> &rhs)
-{
-   Array<Vector *> nullarray(fes.Size());
-   nullarray = NULL;
-
-   ParametricBNLForm::SetParamEssentialBC(bdr_attr_is_ess, nullarray);
-
-   for (int s = 0; s < paramfes.Size(); ++s)
-   {
-      if (rhs[s])
-      {
-         rhs[s]->SetSubVector(*paramess_tdofs[s], 0.0);
-      }
-   }
-}
-
-real_t ParParametricBNLForm::GetEnergy(const Vector &x) const
-{
-   xs_true.Update(const_cast<Vector&>(x), 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));
-   }
-
-   real_t enloc = ParametricBNLForm::GetEnergyBlocked(xs,xdv);
-   real_t englo = 0.0;
-
-   MPI_Allreduce(&enloc, &englo, 1, MPITypeMap<real_t>::mpi_type, MPI_SUM,
-                 ParFESpace(0)->GetComm());
-
-   return englo;
-}
-
-void ParParametricBNLForm::Mult(const Vector &x, Vector &y) const
-{
-   xs_true.Update(const_cast<Vector&>(x), block_trueOffsets);
-   ys_true.Update(y, 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));
-   }
-
-   ParametricBNLForm::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 ParParametricBNLForm::ParamMult(const Vector &x, Vector &y) const
-{
-   xs_true.Update(const_cast<Vector&>(x), paramblock_trueOffsets);
-   ys_true.Update(y, paramblock_trueOffsets);
-   prmxs.Update(paramblock_offsets);
-   prmys.Update(paramblock_offsets);
-
-   for (int s=0; s<paramfes.Size(); ++s)
-   {
-      paramfes[s]->GetProlongationMatrix()->Mult(
-         xs_true.GetBlock(s), prmxs.GetBlock(s));
-   }
-
-   ParametricBNLForm::MultParamBlocked(xsv,adv,xdv,prmys);
-
-   if (fnfi.Size() > 0)
-   {
-      MFEM_ABORT("TODO: assemble contributions from shared face terms");
-   }
-
-   for (int s=0; s<paramfes.Size(); ++s)
-   {
-      paramfes[s]->GetProlongationMatrix()->MultTranspose(
-         prmys.GetBlock(s), ys_true.GetBlock(s));
-
-      ys_true.GetBlock(s).SetSubVector(*paramess_tdofs[s], 0.0);
-   }
-
-}
-
-/// Return the local gradient matrix for the given true-dof vector x
-const BlockOperator & ParParametricBNLForm::GetLocalGradient(
-   const Vector &x) const
-{
-   xs_true.Update(const_cast<Vector&>(x), 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));
-   }
-
-   ParametricBNLForm::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 ParParametricBNLForm::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 & ParParametricBNLForm::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;
-}
-
-ParParametricBNLForm::~ParParametricBNLForm()
-{
-   delete pBlockGrad;
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         delete phBlockGrad(s1,s2);
-      }
-   }
-}
-
-
-void ParParametricBNLForm::SetStateFields(const Vector &xv) const
-{
-   xs_true.Update(const_cast<Vector&>(xv), 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 ParParametricBNLForm::SetAdjointFields(const Vector &av) const
-{
-   xs_true.Update(const_cast<Vector&>(av), 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 ParParametricBNLForm::SetParamFields(const Vector &dv) const
-{
-   xs_true.Update(const_cast<Vector&>(dv),paramblock_trueOffsets);
-   xdv.Update(paramblock_offsets);
-   for (int s=0; s<paramfes.Size(); ++s)
-   {
-      paramfes[s]->GetProlongationMatrix()->Mult(
-         xs_true.GetBlock(s), xdv.GetBlock(s));
-   }
-}
-
-}
-
-#endif

miniapps/mtop/pparamnonlinearform.hpp

@@ -1,114 +0,0 @@
-// Copyright (c) 2010-2025, 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
-
-
-#ifdef MFEM_USE_MPI
-
-#include "mfem.hpp"
-#include "paramnonlinearform.hpp"
-
-namespace mfem
-{
-
-/** @brief A class representing a general parametric parallel block nonlinear
-    operator defined on the Cartesian product of multiple
-    ParFiniteElementSpace%s. */
-/** The ParParametricBNLForm takes as input, and returns as output, vectors on
-    the true dofs. */
-class ParParametricBNLForm : public ParametricBNLForm
-{
-protected:
-   mutable BlockVector xs_true, ys_true;
-   mutable Array2D<OperatorHandle *> phBlockGrad;
-   mutable BlockOperator *pBlockGrad;
-
-public:
-   /// Computes the energy of the system
-   real_t GetEnergy(const Vector &x) const override;
-
-   /// Construct an empty ParParametricBNLForm. Initialize with SetParSpaces().
-   ParParametricBNLForm() : pBlockGrad(nullptr) { }
-
-   /** @brief Construct a ParParametricBNLForm on the given set of
-       parametric and state ParFiniteElementSpace%s. */
-   ParParametricBNLForm(Array<ParFiniteElementSpace *> &statef,
-                        Array<ParFiniteElementSpace *> &paramf);
-
-   /// Return the @a k-th parallel FE state space of the ParParametricBNLForm.
-   ParFiniteElementSpace *ParFESpace(int k);
-   /** @brief Return the @a k-th parallel FE state space of the
-       ParParametricBNLForm (const version). */
-   const ParFiniteElementSpace *ParFESpace(int k) const;
-
-   /// Return the @a k-th parallel FE parameters space of the
-   /// ParParametricBNLForm.
-   ParFiniteElementSpace *ParParamFESpace(int k);
-   /** @brief Return the @a k-th parallel FE parameters space of the
-       ParParametricBNLForm (const version). */
-   const ParFiniteElementSpace *ParParamFESpace(int k) const;
-
-   /** @brief Set the parallel FE spaces for the state and the parametric
-    *  fields. 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 *> &statef,
-                     Array<ParFiniteElementSpace *> &paramf);
-
-   /// Set the state essential BCs. Here, rhs is a true dof vector!
-   void SetEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
-                       Array<Vector *> &rhs) override;
-
-   // Set the essential BCs for the parametric fields. Here, rhs is a true dof
-   // vector!
-   void SetParamEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
-                            Array<Vector *> &rhs) override;
-
-
-   /** @brief Calculates the residual for a state input given by block T-Vector.
-    *  The result is Block T-Vector! The parametric fields should be set in
-    *  advance by calling SetParamFields(). */
-   void Mult(const Vector &x, Vector &y) const override;
-
-   /** @brief Calculates the product of the adjoint field and the derivative of
-    *  the state residual with respect to the parametric fields. The adjoint and
-    *  the state fields should be set in advance by calling SetAdjointFields()
-    *  and SetStateFields(). The input and the result are block T-Vectors!*/
-   void ParamMult(const Vector &x, Vector &y) const override;
-
-   /// Return the local block gradient matrix for the given true-dof vector x
-   const BlockOperator &GetLocalGradient(const Vector &x) const;
-
-   /// Return the block gradient matrix for the given true-dof vector x
-   BlockOperator &GetGradient(const Vector &x) const override;
-
-   /** @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 ~ParParametricBNLForm();
-
-   /// Set the state fields
-   void SetStateFields(const Vector &xv) const override;
-
-   /// Set the adjoint fields
-   void SetAdjointFields(const Vector &av) const override;
-
-   /// Set the parameters/design fields
-   void SetParamFields(const Vector &dv) const override;
-};
-
-}
-
-#endif
-#endif

miniapps/mtop/rand_eigensolver.cpp

@@ -0,0 +1,375 @@
+#include "rand_eigensolver.hpp"
+
+namespace mfem{
+
+void RandomizedSubspaceIteration::Solve()
+{
+
+    if(A==nullptr){return;}
+
+    int myrank=0;
+#ifdef MFEM_USE_MPI
+    MPI_Comm_rank(comm,&myrank);
+#endif
+    std::random_device rd;
+    std::mt19937 generator(rd());
+
+
+    // Create a normal distribution object
+    std::normal_distribution<real_t> distribution(0.0, 1.0);
+
+
+
+    std::vector<Vector> omega; omega.resize(num_modes);
+
+    //populate omega with standard Gaussian RV
+    for(int i=0;i<num_modes;i++){
+        omega[i].SetSize(A->NumCols());
+        Vector& cv=omega[i];
+
+        for(int j=0;j<A->NumCols();j++)
+        {
+            cv[j]=distribution(generator);
+        }
+    }
+
+    if(ess_tdofs!=nullptr){
+        for(int j=0;j<num_modes;j++){
+            Vector& cv=omega[j];
+            for(int i=0;i<ess_tdofs->Size();i++){
+                cv[(*ess_tdofs)[i]]=0.0;
+            }
+        }
+    }
+
+
+    //initialize the modes
+    for(int i=0;i<num_modes;i++){
+        A->Mult(omega[i],modes[i]);
+    }
+
+    real_t gp;
+
+    //Orthogonalize the modes
+    for(int i=0;i<num_modes;i++){
+       for(int j=0;j<i;j++){
+#ifdef MFEM_USE_MPI
+           gp=InnerProduct (comm, modes[i], modes[j]);
+#else
+           gp=InnerProduct (modes[i], modes[j]);
+#endif
+           modes[i].Add(-gp,modes[j]);
+       }
+
+#ifdef MFEM_USE_MPI
+       gp=InnerProduct (comm, modes[i], modes[i]);
+#else
+       gp=InnerProduct (modes[i], modes[i]);
+#endif
+       if(fabs(gp)>std::numeric_limits<real_t>::epsilon()){
+           modes[i]/=sqrt(gp); //scale the vector
+       }
+    }
+
+
+    for(int it=0;it<iter;it++){
+        for(int i=0;i<num_modes;i++){
+            if(symmetric){
+                A->Mult(modes[i],omega[i]);
+            }else{
+                A->MultTranspose(modes[i],omega[i]);
+            }
+
+            for(int j=0;j<i;j++){
+     #ifdef MFEM_USE_MPI
+                gp=InnerProduct (comm, omega[i], omega[j]);
+     #else
+                gp=InnerProduct (omega[i], omega[j]);
+     #endif
+                omega[i].Add(-gp,omega[j]);
+            }
+
+     #ifdef MFEM_USE_MPI
+            gp=InnerProduct (comm, omega[i], omega[i]);
+     #else
+            gp=InnerProduct (omega[i], omega[i]);
+     #endif
+            if(fabs(gp)>std::numeric_limits<real_t>::epsilon()){
+                omega[i]/=sqrt(gp); //scale the vector
+            }
+        }
+
+        //test the product
+        /*
+         if(myrank==0){std::cout<<std::endl;}
+        for(int i=0;i<num_modes;i++){
+            for(int j=0;j<num_modes;j++){
+                real_t bb=InnerProduct(comm,omega[i],modes[j]);
+                if(myrank==0){std::cout<<bb<<" ";}
+
+            }
+            if(myrank==0){std::cout<<std::endl;}
+        }
+        */
+
+
+        for(int i=0;i<num_modes;i++){
+            A->Mult(omega[i],modes[i]);
+            for(int j=0;j<i;j++){
+     #ifdef MFEM_USE_MPI
+                gp=InnerProduct (comm, modes[i], modes[j]);
+     #else
+                gp=InnerProduct (modes[i], modes[j]);
+     #endif
+                modes[i].Add(-gp,modes[j]);
+            }
+
+     #ifdef MFEM_USE_MPI
+            gp=InnerProduct (comm, modes[i], modes[i]);
+     #else
+            gp=InnerProduct (modes[i], modes[i]);
+     #endif
+            if(fabs(gp)>std::numeric_limits<real_t>::epsilon()){
+                modes[i]/=sqrt(gp); //scale the vector
+            }
+        }
+    }
+
+
+}
+
+
+/// Solves the generalized eigenproble using fixed number of iterations
+/// Does not use adaptive check of the precision just run the
+/// specified number of iterations
+void AdaptiveRandomizedGenEig::SolveNA()
+{
+    if(A==nullptr){return;}
+
+    int myrank=0;
+    MPI_Comm_rank(comm,&myrank);
+
+    std::random_device rd;
+    std::mt19937 generator(rd());
+
+    // Create a normal distribution object
+    std::normal_distribution<real_t> distribution(0.0, 1.0);
+
+    std::vector<Vector> omega; omega.resize(num_modes);
+
+    //populate omega with standard Gaussian RV
+    for(int i=0;i<num_modes;i++){
+        omega[i].SetSize(A->NumCols());
+        Vector& cv=omega[i];
+
+        for(int j=0;j<A->NumCols();j++)
+        {
+            cv[j]=distribution(generator);
+        }
+    }
+
+    //initialize the modes
+    for(int i=0;i<num_modes;i++){
+        A->Mult(omega[i],modes[i]);
+        iB->Mult(modes[i], omega[i]);
+    }
+    //Orthogonalize the modes using A inner product
+    //Ortho(A,omega,modes);
+    if(B==nullptr){
+        Ortho(omega,modes);
+    }else{
+        DenseMatrix AA(num_modes);
+        DenseMatrix BB(num_modes);
+
+        for(int i=0;i<num_modes;i++){
+            A->Mult(omega[i],modes[i]);
+            AA(i,i)=InnerProduct(comm, omega[i], modes[i]);
+            for(int j=0;j<i;j++){
+                AA(i,j)=InnerProduct(comm, omega[i], modes[j]);
+                AA(j,i)=AA(i,j);
+            }
+        }
+
+        for(int i=0;i<num_modes;i++){
+            B->Mult(omega[i],modes[i]);
+            BB(i,i)=InnerProduct(comm, omega[i], modes[i]);
+            for(int j=0;j<i;j++){
+                BB(i,j)=InnerProduct(comm, omega[i], modes[j]);
+                BB(j,i)=BB(i,j);
+            }
+        }
+
+        DenseMatrixGeneralizedEigensystem eig(AA,BB,false,true);
+        eig.Eval();
+        DenseMatrix& evecs=eig.RightEigenvectors();
+
+        for(int i=0;i<num_modes;i++){
+            modes[i]=0.0;
+            for(int j=0;j<num_modes;j++){
+                modes[i].Add(evecs(j,i),omega[j]);
+            }
+        }
+    }
+
+
+
+
+    for(int it=0;it<max_iter;it++){
+        for(int i=0;i<num_modes;i++){
+            iB->Mult(modes[i],omega[i]);
+            A->Mult(omega[i], modes[i]);
+        }
+
+        //Ortho(A,modes,omega);
+        Ortho(modes,omega);
+
+        for(int i=0;i<num_modes;i++){
+            A->Mult(omega[i],modes[i]);
+            iB->Mult(modes[i], omega[i]);
+        }
+        //Ortho(A,omega,modes);
+        if(B==nullptr){
+            Ortho(omega,modes);
+        }else{
+            DenseMatrix AA(num_modes);
+            DenseMatrix BB(num_modes);
+
+            for(int i=0;i<num_modes;i++){
+                A->Mult(omega[i],modes[i]);
+                AA(i,i)=InnerProduct(comm, omega[i], modes[i]);
+                for(int j=0;j<i;j++){
+                    AA(i,j)=InnerProduct(comm, omega[i], modes[j]);
+                    AA(j,i)=AA(i,j);
+                }
+            }
+
+            for(int i=0;i<num_modes;i++){
+                B->Mult(omega[i],modes[i]);
+                BB(i,i)=InnerProduct(comm, omega[i], modes[i]);
+                for(int j=0;j<i;j++){
+                    BB(i,j)=InnerProduct(comm, omega[i], modes[j]);
+                    BB(j,i)=BB(i,j);
+                }
+            }
+
+            DenseMatrixGeneralizedEigensystem eig(AA,BB,false,true);
+            eig.Eval();
+            DenseMatrix& evecs=eig.RightEigenvectors();
+            Vector& evals=eig.EigenvaluesRealPart();
+
+            for(int i=0;i<num_modes;i++){
+                modes[i]=0.0;
+                for(int j=0;j<num_modes;j++){
+                    modes[i].Add(evecs(j,i),omega[j]);
+                }
+            }
+        }
+
+    }
+
+    Ortho(A,modes);
+
+
+}
+
+void AdaptiveRandomizedGenEig::Ortho(const Operator* C,
+                                      std::vector<Vector>& vecs)
+{
+    Vector tv(vecs[0]);
+    DenseMatrix AA(vecs.size(),vecs.size());
+    real_t gp;
+    for(size_t i=0;i<vecs.size();i++){
+        C->Mult(vecs[i],tv);
+        for(size_t j=0;j<i;j++){
+            gp=InnerProduct (comm, tv, vecs[j]);
+            AA(i,j)=gp;
+            AA(j,i)=gp;
+        }
+        gp=InnerProduct (comm, tv, vecs[i]);
+        AA(i,i)=gp;
+    }
+
+    Vector sv(vecs.size());
+
+    CholeskyFactors chol(AA.GetData());
+    chol.Factor(AA.NumRows());
+
+    int myrank;
+    MPI_Comm_rank(comm,&myrank);
+
+    for(int i=(vecs.size()-1);i>=0;i=i-1){
+        sv=0.0; sv[i]=1.0;
+        chol.USolve(vecs.size(),1,sv.GetData());
+        (vecs[i])*=sv[i];
+        for(int j=(i-1);j>=0;j--){
+            vecs[i].Add(sv[j],vecs[j]);
+        }
+    }
+}
+
+
+//return the othogonalized vectors in orth
+void AdaptiveRandomizedGenEig::Ortho(const Operator* B,
+                                      std::vector<Vector>& vecs,
+                                      std::vector<Vector>& orth)
+{
+    real_t gp;
+
+
+    DenseMatrix A(vecs.size(),vecs.size());
+    for(size_t i=0;i<vecs.size();i++){
+        orth[i].SetSize(vecs[i].Size());
+        B->Mult(vecs[i],orth[i]);
+        for(size_t j=0;j<i;j++){
+            gp=InnerProduct (comm, orth[i], vecs[j]);
+            A(i,j)=gp;
+            A(j,i)=gp;
+        }
+        gp=InnerProduct (comm, orth[i], vecs[i]);
+        A(i,i)=gp;
+    }
+
+    Vector eval; eval.SetSize(vecs.size());
+    DenseMatrix evec; evec.SetSize(vecs.size(),vecs.size());
+    A.Eigensystem(eval,evec);
+
+    int nn=vecs.size()-1;
+
+    for(size_t i=0;i<vecs.size();i++){
+        orth[nn-i]=real_t(0.0);
+        for(size_t j=0;j<vecs.size();j++){
+            orth[nn-i].Add(evec(j,i)/sqrt(eval(i)),vecs[j]);
+        }
+    }
+}
+
+void AdaptiveRandomizedGenEig::Ortho( std::vector<Vector>& vecs,
+                                      std::vector<Vector>& orth)
+{
+    real_t gp;
+    DenseMatrix A(vecs.size(),vecs.size());
+    for(size_t i=0;i<vecs.size();i++){
+        for(size_t j=0;j<i;j++){
+            gp=InnerProduct (comm, vecs[i], vecs[j]);
+            A(i,j)=gp;
+            A(j,i)=gp;
+        }
+        gp=InnerProduct (comm, vecs[i], vecs[i]);
+        A(i,i)=gp;
+    }
+
+    Vector eval; eval.SetSize(vecs.size());
+    DenseMatrix evec; evec.SetSize(vecs.size(),vecs.size());
+    A.Eigensystem(eval,evec);
+
+    int nn=vecs.size()-1;
+
+    for(size_t i=0;i<vecs.size();i++){
+        orth[nn-i]=real_t(0.0);
+        for(size_t j=0;j<vecs.size();j++){
+            orth[nn-i].Add(evec(j,i)/sqrt(eval(i)),vecs[j]);
+        }
+    }
+}
+
+}

miniapps/mtop/rand_eigensolver.hpp

@@ -0,0 +1,217 @@
+#ifndef RAND_EIGENSOLVER_HPP
+#define RAND_EIGENSOLVER_HPP
+
+#include "mfem.hpp"
+
+
+namespace mfem{
+
+class RandomizedSubspaceIteration
+{
+public:
+
+    RandomizedSubspaceIteration()
+    {
+#ifdef MFEM_USE_MPI
+        comm=MPI_COMM_WORLD;
+#endif
+        num_modes=1;
+        A=nullptr;
+        modes.resize(num_modes);
+        symmetric=false;
+        ess_tdofs=nullptr;
+        iter=1;
+    }
+
+#ifdef MFEM_USE_MPI
+    RandomizedSubspaceIteration(MPI_Comm comm_, bool symm=false)
+    {
+        comm=comm_;
+        num_modes=1;
+        A=nullptr;
+        modes.resize(num_modes);
+        symmetric=symm;
+        ess_tdofs=nullptr;
+        iter=1;
+    }
+#endif
+
+    void SetConstrDOFs(mfem::Array<int>& ess_tdofs_)
+    {
+        ess_tdofs=&ess_tdofs_;
+    }
+
+    void SetNumModes(int num_)
+    {
+        num_modes=num_;
+        modes.resize(num_modes);
+        if(A!=nullptr)
+        {
+            for(int i=0;i<num_modes;i++){
+                modes[i].SetSize(A->NumRows());
+            }
+        }
+    }
+
+    int GetNumModes() const
+    {
+        return num_modes;
+    }
+
+    void SetNumIter(int it_){
+        iter=it_;
+    }
+
+    int GetNumIter() const
+    {
+        return iter;
+    }
+
+    void SetOperator(const Operator& A_, bool symm=false)
+    {
+        A=&A_;
+        for(int i=0;i<num_modes;i++){
+            modes[i].SetSize(A->NumRows());
+        }
+        symmetric=symm;
+    }
+
+    void Solve();
+
+    void GetMode(int i, Vector& q){
+        if((i<num_modes)&&(A!=nullptr))
+        {
+            q=modes[i];
+        }
+    }
+
+    const std::vector<Vector>& GetModes() const
+    {
+        return modes;
+    }
+
+private:
+#ifdef MFEM_USE_MPI
+    MPI_Comm comm;
+#endif
+
+    mfem::Array<int>* ess_tdofs;
+
+    int num_modes;
+    int iter;
+    const Operator* A;
+    std::vector<Vector> modes;
+
+    bool symmetric;
+
+};
+
+class AdaptiveRandomizedGenEig
+{
+public:
+    AdaptiveRandomizedGenEig(MPI_Comm comm_):comm(comm_)
+    {
+        num_modes=1;
+        max_iter=10;
+        eps=1e-6;
+
+        modes.resize(num_modes);
+        B=nullptr;
+        A=nullptr;
+        iB=nullptr;
+    }
+
+    void SetOperators(const Operator& A_, const Operator& iB_)
+    {
+        A=&A_;
+        iB=&iB_;
+
+        for(int i=0;i<num_modes;i++){
+            modes[i].SetSize(A->NumRows());
+        }
+    }
+
+    void SetOperators(const Operator& A_, const Operator& B_, const Operator& iB_)
+    {
+        A=&A_;
+        iB=&iB_;
+        B=&B_;
+
+        for(int i=0;i<num_modes;i++){
+            modes[i].SetSize(A->NumRows());
+        }
+    }
+
+    void SetNumModes(int num_)
+    {
+        num_modes=num_;
+        modes.resize(num_modes);
+        if(A!=nullptr)
+        {
+            for(int i=0;i<num_modes;i++){
+                modes[i].SetSize(A->NumRows());
+            }
+        }
+    }
+
+    int GetNumModes() const
+    {
+        return num_modes;
+    }
+
+    void SetNumIter(int it_){
+        max_iter=it_;
+    }
+
+    int GetNumIter() const
+    {
+        return max_iter;
+    }
+
+    void Solve(bool flag_adaptive=false);
+
+
+    void GetMode(int i, Vector& q){
+        if((i<num_modes)&&(A!=nullptr))
+        {
+            q=modes[i];
+        }
+    }
+
+    const std::vector<Vector>& GetModes() const
+    {
+        return modes;
+    }
+
+    void Ortho(const Operator* C,
+                std::vector<Vector>& vecs);
+
+    void Ortho(const Operator* C,
+                std::vector<Vector>& vecs,
+                std::vector<Vector>& orth);
+
+    void Ortho( std::vector<Vector>& vecs,
+                std::vector<Vector>& orth);
+
+
+    void SolveNA(); //non-addaptive solve
+private:
+    MPI_Comm comm;
+
+    const Operator* A;
+    const Operator* iB;
+    const Operator* B;
+
+    int num_modes;
+    int max_iter;
+    real_t eps;
+
+    std::vector<Vector> modes;
+    Vector evals;
+
+
+};
+
+
+}
+#endif // RAND_EIGENSOLVER_HPP

miniapps/mtop/seqheat.cpp

@@ -1,308 +0,0 @@
-// Copyright (c) 2010-2025, 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.
-//
-//     ----------------------------------------------------------------
-//     SeqHeat Miniapp: Gradients of PDE constrained objective function
-//     ----------------------------------------------------------------
-//                          (Sequential Version)
-//
-// The following example computes the gradients of a specified objective
-// function with respect to parametric fields. The objective function is having
-// the following form f(u(\rho)) where u(\rho) is a solution of a specific state
-// problem (in the example that is the diffusion equation), and \rho is a
-// parametric field discretized by finite elements. The parametric field (also
-// called density in topology optimization) controls the coefficients of the
-// state equation. For the considered case, the density controls the diffusion
-// coefficient within the computational domain.
-//
-// For more information, the users are referred to:
-//
-//    Hinze, M.; Pinnau, R.; Ulbrich, M. & Ulbrich, S.
-//    Optimization with PDE Constraints
-//    Springer Netherlands, 2009
-//
-//    Bendsøe, M. P. & Sigmund, O.
-//    Topology Optimization - Theory, Methods and Applications
-//    Springer Verlag, Berlin Heidelberg, 2003
-//
-// Compile with: make seqheat
-//
-// Sample runs:
-//
-//    seqheat -m ../../data/star-mixed.mesh
-//    seqheat --visualization
-
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-
-#include "mtop_integrators.hpp"
-
-using namespace mfem;
-
-int main(int argc, char *argv[])
-{
-   const char *mesh_file = "../../data/star.vtk";
-   int ser_ref_levels = 1;
-   int order = 2;
-   bool visualization = false;
-   real_t newton_rel_tol = 1e-4;
-   real_t 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);
-
-   // 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();
-
-   // 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();
-   }
-
-   // Diffusion coefficient
-   mfem::ConstantCoefficient* diffco=new mfem::ConstantCoefficient(1.0);
-   // Heat source
-   mfem::ConstantCoefficient* loadco=new mfem::ConstantCoefficient(1.0);
-   // Define the q-function
-   mfem::QLinearDiffusion* qfun=new mfem::QLinearDiffusion(*diffco,*loadco,1.0,
-                                                           1e-7,4.0,0.5);
-
-   // Define FE collection and space for the state solution
-   mfem::H1_FECollection sfec(order, dim);
-   mfem::FiniteElementSpace* sfes=new mfem::FiniteElementSpace(mesh,&sfec,1);
-   // Define FE collection and space for the density field
-   mfem::L2_FECollection pfec(order, dim);
-   mfem::FiniteElementSpace* pfes=new mfem::FiniteElementSpace(mesh,&pfec,1);
-
-   // Define the arrays for the nonlinear form
-   mfem::Array<mfem::FiniteElementSpace*> asfes;
-   mfem::Array<mfem::FiniteElementSpace*> apfes;
-
-   asfes.Append(sfes);
-   apfes.Append(pfes);
-   // Define parametric block nonlinear form using single scalar H1 field
-   // and L2 scalar density field
-   mfem::ParametricBNLForm* nf=new mfem::ParametricBNLForm(asfes,apfes);
-   // Add the parametric integrator
-   nf->AddDomainIntegrator(new mfem::ParametricLinearDiffusion(*qfun));
-
-   // Define true block vectors for state, adjoint, residual
-   mfem::BlockVector solbv; solbv.Update(nf->GetBlockTrueOffsets());    solbv=0.0;
-   mfem::BlockVector adjbv; adjbv.Update(nf->GetBlockTrueOffsets());    adjbv=0.0;
-   mfem::BlockVector resbv; resbv.Update(nf->GetBlockTrueOffsets());    resbv=0.0;
-   // Define true block vectors for parametric field and gradients
-   mfem::BlockVector prmbv; prmbv.Update(nf->ParamGetBlockTrueOffsets());
-   prmbv=0.0;
-   mfem::BlockVector grdbv; grdbv.Update(nf->ParamGetBlockTrueOffsets());
-   grdbv=0.0;
-
-   // 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);
-   delete ess_bdr[0];
-
-   // Define the linear solvers
-   mfem::GMRESSolver *gmres;
-   gmres = new mfem::GMRESSolver();
-   gmres->SetAbsTol(newton_abs_tol/10);
-   gmres->SetRelTol(newton_rel_tol/10);
-   gmres->SetMaxIter(300);
-   gmres->SetPrintLevel(print_level);
-
-   // Define the Newton solver
-   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);
-
-   // Solve the problem
-   // Set the density to 0.5
-   prmbv=0.5;
-   nf->SetParamFields(prmbv); // Set the density
-   // Define the RHS
-   mfem::Vector b;
-   solbv=0.0;
-   // Newton solve
-   ns->Mult(b, solbv);
-
-   // Compute the residual
-   nf->Mult(solbv,resbv);
-   std::cout<<"Norm residual="<<resbv.Norml2()<<std::endl;
-
-   // Compute the energy of the state system
-   real_t energy = nf->GetEnergy(solbv);
-   std::cout<<"energy ="<< energy<<std::endl;
-
-   // Define the block nonlinear form utilized for representing the
-   // objective. The input is the state array asfes defined earlier.
-   mfem::BlockNonlinearForm* ob=new mfem::BlockNonlinearForm(asfes);
-
-   // Add the integrator for the objective
-   ob->AddDomainIntegrator(new mfem::DiffusionObjIntegrator());
-
-   // Compute the objective
-   real_t obj=ob->GetEnergy(solbv);
-   std::cout<<"Objective ="<<obj<<std::endl;
-
-   // Solve the adjoint
-   {
-      mfem::BlockVector adjrhs; adjrhs.Update(nf->GetBlockTrueOffsets());  adjrhs=0.0;
-      // Compute the RHS for the adjoint
-      ob->Mult(solbv, adjrhs);
-      // Get the tangent matrix from the state problem
-      mfem::BlockOperator& A=nf->GetGradient(solbv);
-      // We do not need to transpose the operator for diffusion
-      gmres->SetOperator(A.GetBlock(0,0));
-      // Compute the adjoint solution
-      gmres->Mult(adjrhs.GetBlock(0), adjbv.GetBlock(0));
-   }
-
-   // Compute gradients
-   nf->SetAdjointFields(adjbv);
-   nf->SetStateFields(solbv);
-   nf->ParamMult(prmbv, grdbv);
-
-   // Dump out the data
-   if (visualization)
-   {
-      mfem::ParaViewDataCollection *dacol=new mfem::ParaViewDataCollection("SeqHeat",
-                                                                           mesh);
-      mfem::GridFunction gfgrd(pfes); gfgrd.SetFromTrueDofs(grdbv.GetBlock(0));
-      mfem::GridFunction gfdns(pfes); gfdns.SetFromTrueDofs(prmbv.GetBlock(0));
-      // Define state grid function
-      mfem::GridFunction gfsol(sfes); gfsol.SetFromTrueDofs(solbv.GetBlock(0));
-      mfem::GridFunction gfadj(sfes); gfadj.SetFromTrueDofs(adjbv.GetBlock(0));
-
-      dacol->SetLevelsOfDetail(order);
-      dacol->RegisterField("sol", &gfsol);
-      dacol->RegisterField("adj", &gfadj);
-      dacol->RegisterField("dns", &gfdns);
-      dacol->RegisterField("grd", &gfgrd);
-
-      dacol->SetTime(1.0);
-      dacol->SetCycle(1);
-      dacol->Save();
-
-      delete dacol;
-   }
-
-   // FD check
-   {
-      // Perturbation vector
-      mfem::BlockVector prtbv;
-      mfem::BlockVector tmpbv;
-      prtbv.Update(nf->ParamGetBlockTrueOffsets());
-      tmpbv.Update(nf->ParamGetBlockTrueOffsets());
-      // Generate the perturbation
-      prtbv.GetBlock(0).Randomize();
-      prtbv*=1.0;
-      // Scaling parameter
-      real_t lsc=1.0;
-
-      // Compute initial objective
-      real_t gQoI=ob->GetEnergy(solbv);
-      real_t lQoI;
-
-      // Norm of the perturbation
-      real_t nd=mfem::InnerProduct(prtbv,prtbv);
-      // Projection of the adjoint gradient on the perturbation
-      real_t td=mfem::InnerProduct(prtbv,grdbv);
-      // Normalize the directional derivative
-      td=td/nd;
-
-      for (int l = 0; l < 10; l++)
-      {
-         lsc/=10.0;
-         // Scale the perturbation
-         prtbv/=10.0;
-         // Add the perturbation to the original density
-         add(prmbv,prtbv,tmpbv);
-         nf->SetParamFields(tmpbv);
-         // Solve the physics
-         ns->Mult(b,solbv);
-         // Compute the objective
-         lQoI=ob->GetEnergy(solbv);
-         // FD approximation
-         real_t ld=(lQoI-gQoI)/lsc;
-         std::cout << "dx=" << lsc << " FD gradient=" << ld/nd
-                   << " adjoint gradient=" << td
-                   << " err=" << std::fabs(ld/nd-td) << std::endl;
-      }
-   }
-
-   delete ob;
-
-   delete ns;
-   delete gmres;
-
-   delete nf;
-   delete pfes;
-   delete sfes;
-
-   delete qfun;
-   delete loadco;
-   delete diffco;
-
-   delete mesh;
-
-   return 0;
-}

miniapps/mtop/stokes_solver.cpp

@@ -0,0 +1,579 @@
+#include "stokes_solver.hpp"
+#include "general/forall.hpp"
+
+StokesOperator::StokesOperator(mfem::ParMesh* mesh_,int vorder):zeroc(0.0),onec(1.0)
+{
+    if(vorder<2){vorder=2;}
+
+    pmesh=mesh_;
+    const int dim=pmesh->Dimension();
+    vfec=new mfem::H1_FECollection(vorder,dim);
+    vfes=new mfem::ParFiniteElementSpace(pmesh,vfec,dim, mfem::Ordering::byVDIM);
+
+    pfec=new mfem::H1_FECollection(vorder-1,dim);
+    pfes=new mfem::ParFiniteElementSpace(pmesh,pfec,1, mfem::Ordering::byVDIM);
+
+    fvelo.SetSpace(vfes); fvelo=0.0;
+    avelo.SetSpace(vfes); avelo=0.0;
+
+    fpres.SetSpace(pfes); fpres=0.0;
+    apres.SetSpace(pfes); apres=0.0;
+
+    block_offsets.SetSize(3);
+    block_offsets[0] = 0;
+    block_offsets[1] = vfes->GetVSize();
+    block_offsets[2] = pfes->GetVSize();
+    block_offsets.PartialSum();
+
+    block_true_offsets.SetSize(3);
+    block_true_offsets[0] = 0;
+    block_true_offsets[1] = vfes->TrueVSize();
+    block_true_offsets[2] = pfes->TrueVSize();
+    block_true_offsets.PartialSum();
+
+    sol.Update(block_true_offsets);
+    adj.Update(block_true_offsets);
+    rhs.Update(block_true_offsets);
+
+    SetLinearSolver();
+
+    mfem::Operator::width=block_true_offsets[2];
+    mfem::Operator::height=block_true_offsets[2];
+
+    mu=nullptr;
+    brink=nullptr;
+
+    af=nullptr;
+    mf=nullptr;
+
+    lf=nullptr;
+    pf=nullptr;
+
+
+    bf=new mfem::ParMixedBilinearForm(vfes,pfes);
+    bf->AddDomainIntegrator(new mfem::VectorDivergenceIntegrator());
+    bf->Assemble(0);
+    bf->Finalize();
+
+
+    bop=new mfem::BlockOperator(block_true_offsets); bop->owns_blocks=false;
+
+    pop=nullptr; //preconditioner
+    ls=nullptr;  //linear solver
+}
+
+
+StokesOperator::~StokesOperator()
+{
+    delete bop;
+
+    delete pop;
+    delete ls;
+
+    delete pf;
+    delete lf;
+    delete bf;
+    delete af;
+
+    delete pfes;
+    delete pfec;
+    delete vfes;
+    delete vfec;
+}
+
+void StokesOperator::SetLinearSolver(mfem::real_t rtol, mfem::real_t atol, int miter)
+{
+    linear_rtol=rtol;
+    linear_atol=atol;
+    linear_iter=miter;
+}
+
+void StokesOperator::AddVelocityBC(int id, int dir, mfem::real_t val)
+{
+    if(dir==0){
+        bcx[id]=mfem::ConstantCoefficient(val);
+        AddVelocityBC(id,dir,bcx[id]);
+    }
+    if(dir==1){
+        bcy[id]=mfem::ConstantCoefficient(val);
+        AddVelocityBC(id,dir,bcy[id]);
+
+    }
+    if(dir==2){
+        bcz[id]=mfem::ConstantCoefficient(val);
+        AddVelocityBC(id,dir,bcz[id]);
+    }
+    if(dir==4){
+        bcx[id]=mfem::ConstantCoefficient(val);
+        bcy[id]=mfem::ConstantCoefficient(val);
+        bcz[id]=mfem::ConstantCoefficient(val);
+        AddVelocityBC(id,0,bcx[id]);
+        AddVelocityBC(id,1,bcy[id]);
+        AddVelocityBC(id,2,bcz[id]);
+    }
+}
+
+void StokesOperator::DelVelocityBC()
+{
+    bccx.clear();
+    bccy.clear();
+    bccz.clear();
+
+    bcx.clear();
+    bcy.clear();
+    bcz.clear();
+
+    ess_tdofv.DeleteAll();
+}
+
+void StokesOperator::AddVelocityBC(int id, int dir, mfem::Coefficient &val)
+{
+    if(dir==0){ bccx[id]=&val; }
+    if(dir==1){ bccy[id]=&val; }
+    if(dir==2){ bccz[id]=&val; }
+    if(dir==4){ bccx[id]=&val; bccy[id]=&val; bccz[id]=&val;}
+    if(pmesh->Dimension()==2)
+    {
+        bccz.clear();
+    }
+}
+
+void StokesOperator::SetEssTDofs(mfem::Vector& bsol, mfem::Array<int>& ess_dofs)
+{
+    // Set the BC
+
+    ess_tdofv.DeleteAll();
+
+    mfem::Array<int> ess_tdofx;
+    mfem::Array<int> ess_tdofy;
+    mfem::Array<int> ess_tdofz;
+
+    int dim=pmesh->Dimension();
+
+    {
+        for(auto it=bccx.begin();it!=bccx.end();it++)
+        {
+            mfem::Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+            ess_bdr=0;
+            ess_bdr[it->first -1]=1;
+            mfem::Array<int> ess_tdof_list;
+            vfes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list,0);
+            ess_tdofx.Append(ess_tdof_list);
+
+            mfem::VectorArrayCoefficient pcoeff(dim);
+            pcoeff.Set(0, it->second, false);
+            fvelo.ProjectBdrCoefficient(pcoeff, ess_bdr);
+        }
+        //copy tdofsx from velocity grid function
+        {
+            mfem::Vector& vc=fvelo.GetTrueVector();
+            for(int ii=0;ii<ess_tdofx.Size();ii++)
+            {
+                bsol[ess_tdofx[ii]]=vc[ess_tdofx[ii]];
+            }
+        }
+        ess_dofs.Append(ess_tdofx); ess_tdofx.DeleteAll();
+
+        for(auto it=bccy.begin();it!=bccy.end();it++)
+        {
+            mfem::Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+            ess_bdr=0;
+            ess_bdr[it->first -1]=1;
+            mfem::Array<int> ess_tdof_list;
+            vfes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list,1);
+            ess_tdofy.Append(ess_tdof_list);
+
+            mfem::VectorArrayCoefficient pcoeff(dim);
+            pcoeff.Set(1, it->second, false);
+            fvelo.ProjectBdrCoefficient(pcoeff, ess_bdr);
+        }
+        //copy tdofsy from velocity grid function
+        {
+            mfem::Vector& vc=fvelo.GetTrueVector();
+            for(int ii=0;ii<ess_tdofy.Size();ii++)
+            {
+                bsol[ess_tdofy[ii]]=vc[ess_tdofy[ii]];
+            }
+        }
+        ess_dofs.Append(ess_tdofy); ess_tdofy.DeleteAll();
+
+        if(dim==3){
+            for(auto it=bccz.begin();it!=bccz.end();it++)
+            {
+                mfem::Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+                ess_bdr=0;
+                ess_bdr[it->first -1]=1;
+                mfem::Array<int> ess_tdof_list;
+                vfes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list,2);
+                ess_tdofz.Append(ess_tdof_list);
+
+                mfem::VectorArrayCoefficient pcoeff(dim);
+                pcoeff.Set(2, it->second, false);
+                fvelo.ProjectBdrCoefficient(pcoeff, ess_bdr);
+            }
+
+            //copy tdofsz from velocity grid function
+            {
+                mfem::Vector& vc=fvelo.GetTrueVector();
+                for(int ii=0;ii<ess_tdofz.Size();ii++)
+                {
+                    bsol[ess_tdofz[ii]]=vc[ess_tdofz[ii]];
+                }
+            }
+            ess_dofs.Append(ess_tdofz); ess_tdofz.DeleteAll();
+        }
+    }
+}
+
+
+void StokesOperator::AssemblePrec1()
+{
+    delete pop;
+    mfem::BlockDiagonalPreconditioner* dpop=new mfem::BlockDiagonalPreconditioner(block_true_offsets);
+    dpop->owns_blocks=true;
+
+    mfem::HypreBoomerAMG* prec=new mfem::HypreBoomerAMG();
+    prec->SetElasticityOptions(vfes);
+    prec->SetPrintLevel(1);
+    prec->SetOperator(*A);
+
+    //dpop->SetDiagonalBlock(0,new mfem::IdentityOperator(vfes->TrueVSize()));
+    dpop->SetDiagonalBlock(0,prec);
+    dpop->SetDiagonalBlock(1,new mfem::IdentityOperator(pfes->TrueVSize()));
+    pop=dpop;
+}
+
+void StokesOperator::AssemblePrec2()
+{
+    delete pop;
+    mfem::BlockDiagonalPreconditioner* dpop=new mfem::BlockDiagonalPreconditioner(block_true_offsets);
+    dpop->owns_blocks=true;
+
+
+    mfem::HypreBoomerAMG* prec=new mfem::HypreBoomerAMG();
+    prec->SetElasticityOptions(vfes);
+    prec->SetPrintLevel(1);
+    prec->SetOperator(*A);
+
+    ivisc=std::unique_ptr<InverseCoeff>(new InverseCoeff(mu));
+    mf=std::unique_ptr<mfem::ParBilinearForm>(new mfem::ParBilinearForm(pfes));
+    mf->AddDomainIntegrator(new mfem::MassIntegrator(*ivisc));
+    mf->Assemble(0);
+    mf->Finalize();
+
+    M=std::unique_ptr<mfem::HypreParMatrix>(mf->ParallelAssemble());
+    mfem::HypreBoomerAMG* prem=new mfem::HypreBoomerAMG();
+    prem->SetPrintLevel(1);
+    prem->SetOperator(*M);
+
+    dpop->SetDiagonalBlock(0,prec);
+    dpop->SetDiagonalBlock(1,prem);
+
+    pop=dpop;
+}
+
+void StokesOperator::AssemblePrec4()
+{
+    delete pop;
+    mfem::BlockLowerTriangularPreconditioner* dpop=new mfem::BlockLowerTriangularPreconditioner(block_true_offsets);
+    dpop->owns_blocks=false;
+
+    prec1=std::unique_ptr<mfem::HypreBoomerAMG>(new mfem::HypreBoomerAMG());
+    //prec1->SetElasticityOptions(vfes);
+    prec1->SetPrintLevel(1);
+    prec1->SetOperator(*A);
+
+    ivisc=std::unique_ptr<InverseCoeff>(new InverseCoeff(mu));
+
+    mf=std::unique_ptr<mfem::ParBilinearForm>(new mfem::ParBilinearForm(pfes));
+    mf->AddDomainIntegrator(new mfem::MassIntegrator(*ivisc));
+    mf->Assemble(0);
+    mf->Finalize();
+
+    M=std::unique_ptr<mfem::HypreParMatrix>(mf->ParallelAssemble());
+    prec2=std::unique_ptr<mfem::HypreBoomerAMG>(new mfem::HypreBoomerAMG());
+    prec2->SetPrintLevel(1);
+    prec2->SetOperator(*M);
+
+    dpop->SetBlock(1,0,B.get());
+    dpop->SetDiagonalBlock(0,prec1.get());
+    dpop->SetDiagonalBlock(1,prec2.get());
+
+    pop=dpop;
+}
+
+void StokesOperator::AssemblePrec3()
+{
+    delete pop;
+    mfem::BlockDiagonalPreconditioner* dpop=new mfem::BlockDiagonalPreconditioner(block_true_offsets);
+    dpop->owns_blocks=true;
+
+    mfem::CGSolver*  ls1=new mfem::CGSolver(pmesh->GetComm());
+    ls1->SetAbsTol(linear_atol);
+    ls1->SetRelTol(linear_rtol);
+    ls1->SetMaxIter(100);
+
+    prec1=std::unique_ptr<mfem::HypreBoomerAMG>(new mfem::HypreBoomerAMG());
+    //set the rigid body modes
+    prec1->SetElasticityOptions(vfes);
+    prec1->SetPrintLevel(1);
+    ls1->SetPreconditioner(*prec1);
+    ls1->SetOperator(*A);
+    ls1->SetPrintLevel(0);
+
+
+    ivisc=std::unique_ptr<InverseCoeff>(new InverseCoeff(mu));
+
+    mf=std::unique_ptr<mfem::ParBilinearForm>(new mfem::ParBilinearForm(pfes));
+    mf->AddDomainIntegrator(new mfem::MassIntegrator(*ivisc));
+    mf->Assemble(0);
+    mf->Finalize();
+
+    M=std::unique_ptr<mfem::HypreParMatrix>(mf->ParallelAssemble());
+
+    mfem::CGSolver*  ls2=new mfem::CGSolver(pmesh->GetComm());
+    ls2->SetAbsTol(linear_atol);
+    ls2->SetRelTol(linear_rtol);
+    ls2->SetMaxIter(100);
+
+    prec2=std::unique_ptr<mfem::HypreBoomerAMG>(new mfem::HypreBoomerAMG());
+    prec2->SetPrintLevel(1);
+    ls2->SetPreconditioner(*prec2);
+    ls2->SetOperator(*M);
+    ls2->SetPrintLevel(0);
+
+
+    dpop->SetDiagonalBlock(0,ls1);
+    dpop->SetDiagonalBlock(1,ls2);
+
+    pop=dpop;
+
+}
+
+void StokesOperator::AssemblePrec5()
+{
+    //get the diagonal of A
+    mfem::HypreParVector* Ad = new mfem::HypreParVector(pmesh->GetComm(),
+                                                        A->GetGlobalNumRows(),
+                                                        A->GetRowStarts());
+    A->GetDiag(*Ad);
+
+    //build approximation to the Schur complement
+    std::unique_ptr<mfem::HypreParMatrix> MBt(B->Transpose());
+    MBt->InvScaleRows(*Ad);
+    M=std::unique_ptr<mfem::HypreParMatrix>(mfem::ParMult(B.get(),MBt.get()));
+
+    //prepare the block struture
+    delete pop;
+    mfem::BlockDiagonalPreconditioner* dpop=new mfem::BlockDiagonalPreconditioner(block_true_offsets);
+    dpop->owns_blocks=true;
+
+
+    mfem::HypreBoomerAMG* prec=new mfem::HypreBoomerAMG();
+    prec->SetElasticityOptions(vfes);
+    prec->SetPrintLevel(1);
+    prec->SetOperator(*A);
+
+    mfem::HypreBoomerAMG* prem=new mfem::HypreBoomerAMG();
+    prem->SetPrintLevel(1);
+    prem->SetOperator(*M);
+
+    dpop->SetDiagonalBlock(0,prec);
+    dpop->SetDiagonalBlock(1,prem);
+
+    pop=dpop;
+
+    delete Ad;
+
+}
+
+void StokesOperator::AssemblePrec6()
+{
+    M=std::unique_ptr<mfem::HypreParMatrix>(mfem::ParMult(B.get(),D.get()));
+    LSC* prem=new LSC(*B,*A,*D,*M);
+
+    mfem::HypreBoomerAMG* prec=new mfem::HypreBoomerAMG();
+    prec->SetElasticityOptions(vfes);
+    prec->SetPrintLevel(1);
+    prec->SetOperator(*A);
+
+    //prepare the block struture
+    delete pop;
+    mfem::BlockDiagonalPreconditioner* dpop=new mfem::BlockDiagonalPreconditioner(block_true_offsets);
+    dpop->owns_blocks=true;
+
+    dpop->SetDiagonalBlock(0,prec);
+    dpop->SetDiagonalBlock(1,prem);
+
+    pop=dpop;
+}
+
+void StokesOperator::AssemblePrec7()
+{
+    delete pop;
+    mfem::BlockLowerTriangularPreconditioner* dpop=new mfem::BlockLowerTriangularPreconditioner(block_true_offsets);
+    dpop->owns_blocks=false;
+
+    prec1=std::unique_ptr<mfem::HypreBoomerAMG>(new mfem::HypreBoomerAMG());
+    prec1->SetElasticityOptions(vfes);
+    prec1->SetPrintLevel(1);
+    prec1->SetOperator(*A);
+
+    ivisc=std::unique_ptr<InverseCoeff>(new InverseCoeff(mu));
+
+    mf=std::unique_ptr<mfem::ParBilinearForm>(new mfem::ParBilinearForm(pfes));
+    mf->AddDomainIntegrator(new mfem::MassIntegrator(*ivisc));
+    mf->Assemble(0);
+    mf->Finalize();
+
+    M=std::unique_ptr<mfem::HypreParMatrix>(mf->ParallelAssemble());
+    prec2=std::unique_ptr<mfem::HypreBoomerAMG>(new mfem::HypreBoomerAMG());
+    prec2->SetPrintLevel(1);
+    prec2->SetOperator(*M);
+
+    dpop->SetBlock(1,0,B.get());
+    dpop->SetDiagonalBlock(0,prec1.get());
+    dpop->SetDiagonalBlock(1,prec2.get());
+
+    pop=dpop;
+}
+
+
+void StokesOperator::Assemble()
+{
+   if(mu==nullptr){return;}
+   //allocate af
+   if(af==nullptr){
+      af=new mfem::ParBilinearForm(vfes);
+      af->SetDiagonalPolicy(DIAG_ONE);
+      af->AddDomainIntegrator(new mfem::ElasticityIntegrator(zeroc,*mu));
+      if(brink!=nullptr){
+          af->AddDomainIntegrator(new mfem::VectorMassIntegrator(*brink));
+      }
+   }
+
+   //set BC
+   sol=mfem::real_t(0.0);
+   SetEssTDofs(sol.GetBlock(0),ess_tdofv); //set BC for the velocity
+
+   af->Assemble(0);
+   af->Finalize();
+   //af->FormSystemMatrix(ess_tdofv,A);
+   //std::unique_ptr<mfem::HypreParMatrix> Ael(af->ParallelAssembleElim()); Ae=std::move(Ael);
+   A=std::unique_ptr<mfem::HypreParMatrix>(af->ParallelAssemble());
+   Ae=std::unique_ptr<mfem::HypreParMatrix>(A->EliminateRowsCols(ess_tdofv));
+
+   std::cout<<"H="<<A->Height()<<" W="<<A->Width()<<std::endl;
+
+   B=std::unique_ptr<mfem::HypreParMatrix>(bf->ParallelAssemble());
+   std::cout<<"H="<<B->Height()<<" W="<<B->Width()<<std::endl;
+
+   D=std::unique_ptr<mfem::HypreParMatrix>(B->Transpose());
+   std::cout<<"H="<<D->Height()<<" W="<<D->Width()<<std::endl;
+
+   B->EliminateRows(ess_tdofp);
+   Be=std::unique_ptr<mfem::HypreParMatrix>(B->EliminateCols(ess_tdofv));
+   D->EliminateRows(ess_tdofv);
+   De=std::unique_ptr<mfem::HypreParMatrix>(D->EliminateCols(ess_tdofp));
+
+   bop->SetBlock(0,0,A.get());
+   bop->SetBlock(0,1,D.get());
+   bop->SetBlock(1,0,B.get());
+
+   //assemble the preconditioner
+   //AssemblePrec1();
+   //AssemblePrec2();
+   //AssemblePrec3();
+   //AssemblePrec4();
+   //AssemblePrec5();
+   AssemblePrec6();
+
+   if(ls==nullptr){
+       ls=new mfem::MINRESSolver(pmesh->GetComm());
+       //ls=new mfem::GMRESSolver(pmesh->GetComm());
+       ls->SetAbsTol(linear_atol);
+       ls->SetRelTol(linear_rtol);
+       ls->SetMaxIter(linear_iter);
+
+       ls->SetOperator(*bop);
+       ls->SetPreconditioner(*pop);
+       ls->SetPrintLevel(1);
+   }else{
+       ls->SetOperator(*bop);
+       ls->SetPreconditioner(*pop);
+   }
+}
+
+void StokesOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const
+{
+    //the rhs x is assumed to have the contribution of the BC
+    ls->Mult(x,y);
+
+    int N=ess_tdofv.Size();
+    mfem::real_t *yp = y.ReadWrite();
+    const mfem::real_t *sp = sol.Read();
+    const int *ep = ess_tdofv.Read();
+    mfem::forall(N, [=] MFEM_HOST_DEVICE (int i)
+    {
+        yp[ep[i]]=sp[ep[i]];
+    });
+
+    //add the pressure BC
+}
+
+void StokesOperator:: MultTranspose(const mfem::Vector &x, mfem::Vector &y) const
+{
+    //the adjoint rhs is assumed to be corrected for the BC
+    //K is symmetric
+    ls->Mult(x,y);
+
+    int N=ess_tdofv.Size();
+    ess_tdofv.Read();
+
+    mfem::real_t *yp = y.ReadWrite();
+    const int *ep = ess_tdofv.Read();
+
+    mfem::forall(N,[=] MFEM_HOST_DEVICE (int i)
+    {
+        yp[ep[i]]=mfem::real_t(0.0);
+    });
+
+    //add the pressure BC
+}
+
+void StokesOperator::FSolve()
+{
+
+    rhs=mfem::real_t(0.0);
+
+    if(lf==nullptr){
+        lf=new mfem::ParLinearForm(vfes);
+        /*
+        if(volforce!=nullptr){
+            lf->AddDomainIntegrator(new mfem::VectorDomainLFIntegrator(*volforce));
+        }
+        //add surface loads
+        */
+
+    }
+
+    (*lf)=mfem::real_t(0.0);
+
+    lf->Assemble();
+    lf->ParallelAssemble(rhs.GetBlock(0));
+
+
+
+    Ae->AddMult(sol.GetBlock(0),rhs.GetBlock(0),-1.0);
+    Be->AddMult(sol.GetBlock(0),rhs.GetBlock(1),-1.0);
+    De->AddMult(sol.GetBlock(1),rhs.GetBlock(0),-1.0);
+
+    /*
+    for(int i=0;i<ess_tdofv.Size();i++){
+        rhs.GetBlock(0)[ess_tdofv[i]]=sol.GetBlock(0)[ess_tdofv[i]];
+    }
+    */
+
+    ls->Mult(rhs,sol);
+
+}

miniapps/mtop/stokes_solver.hpp

@@ -0,0 +1,307 @@
+#ifndef STOKES_SOLVER_HPP
+#define STOKES_SOLVER_HPP
+
+#include "mfem.hpp"
+
+class StokesOperator:public mfem::Operator
+{
+public:
+    StokesOperator(mfem::ParMesh* mesh_,int vorder=2);
+
+    virtual ~StokesOperator();
+
+    /// Set the Linear Solver
+    void SetLinearSolver(double rtol=1e-8, double atol=1e-12, int miter=1000);
+
+    /// Solves the forward problem.
+    void FSolve();
+
+    /// Forms the tangent matrix
+    void AssembleTangent();
+
+    /// Solves the adjoint with the provided rhs.
+    void ASolve(mfem::Vector& rhs);
+
+    /// Solves the forward problem with the provided rhs.
+    void FSolve(mfem::Vector& rhs);
+
+    /// Adds velocity BC in direction 0(x),1(y),2(z), or 4(all).
+    void AddVelocityBC(int id, int dir, double val);
+
+    /// Adds velocity BC in direction 0(x),1(y),2(z), or 4(all).
+    void AddVelocityBC(int id, int dir, mfem::Coefficient& val);
+
+    /// Clear all velocity BC
+    void DelVelocityBC();
+
+    /// Set the values of the volumetric force.
+    void SetVolForce(double fx,double fy, double fz=0.0);
+
+    /// Add surface load
+    void AddSurfLoad(int id, double fx,double fy, double fz=0.0)
+    {
+        /*
+        mfem::Vector vec; vec.SetSize(pmesh->SpaceDimension());
+        vec[0]=fx;
+        vec[1]=fy;
+        if(pmesh->SpaceDimension()==3){vec[2]=fz;}
+        mfem::VectorConstantCoefficient* vc=new mfem::VectorConstantCoefficient(vec);
+        if(load_coeff.find(id)!=load_coeff.end()){ delete load_coeff[id];}
+        load_coeff[id]=vc;
+        */
+    }
+
+    /// Add surface load
+    void AddSurfLoad(int id, mfem::VectorCoefficient& ff)
+    {
+        /*
+        surf_loads[id]=&ff;
+        */
+    }
+
+    /// Associates coefficient to the volumetric force.
+    void SetVolForce(mfem::VectorCoefficient& ff);
+
+
+    /// Sets BC dofs, bilinear form, preconditioner and solver.
+    /// Should be called before calling Mult of MultTranspose
+    virtual void Assemble();
+
+    /// Forward solve with given RHS. x is the RHS vector. The BC are set to zero.
+    virtual void Mult(const mfem::Vector &x, mfem::Vector &y) const override;
+
+    /// Adjoint solve with given RHS. x is the RHS vector. The BC are set to zero.
+    virtual void MultTranspose(const mfem::Vector &x, mfem::Vector &y) const override;
+
+    /// Set material
+    void SetViscosoty(mfem::Coefficient& mu_)
+    {
+        mu=&mu_;
+
+        delete af; af=nullptr;
+    }
+
+    void SetBrinkman(mfem::Coefficient& alpha_)
+    {
+        brink=&alpha_;
+
+        delete af; af=nullptr;
+    }
+
+    /// Returns the velocity field
+    mfem::ParGridFunction& GetVelocity()
+    {
+        fvelo.SetFromTrueDofs(sol.GetBlock(0));
+        return fvelo;
+    }
+
+    /// Return the pressure field
+    mfem::ParGridFunction& GetPressure()
+    {
+        fpres.SetFromTrueDofs(sol.GetBlock(1));
+        return fpres;
+    }
+
+    void GetVelocity(mfem::ParGridFunction& v)
+    {
+        v.SetSpace(vfes);
+        v.SetFromTrueDofs(sol.GetBlock(0));
+    }
+
+    void GetPressure(mfem::ParGridFunction& p)
+    {
+        p.SetSpace(pfes);
+        p.SetFromTrueDofs(sol.GetBlock(1));
+    }
+
+
+    class InverseCoeff:public mfem::Coefficient
+    {
+    public:
+        InverseCoeff(mfem::Coefficient* co)
+        {
+            cc=co;
+        }
+
+        void SetCoeff(mfem::Coefficient* co)
+        {
+            cc=co;
+        }
+
+
+        virtual mfem::real_t
+        Eval(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip) override
+        {
+            mfem::real_t vv=cc->Eval(T,ip);
+            return 2.0/(1.0*vv);
+            //return vv;
+        }
+
+    private:
+        mfem::Coefficient* cc;
+    };
+
+protected:
+    mfem::ParMesh* pmesh;
+
+    mfem::Array<int> block_offsets; // number of variables + 1
+    mfem::Array<int> block_true_offsets;
+
+    //block vectors
+    mutable mfem::BlockVector sol;
+    mutable mfem::BlockVector adj;
+    mutable mfem::BlockVector rhs;
+
+    mfem::ParGridFunction fvelo;
+    mfem::ParGridFunction fpres;
+
+    mfem::ParGridFunction avelo;
+    mfem::ParGridFunction apres;
+
+    //Linear solver parameters
+    double linear_rtol;
+    double linear_atol;
+    int linear_iter;
+
+    //finite element space for velocity
+    mfem::ParFiniteElementSpace* vfes;
+    //finite element collection for velocity
+    mfem::FiniteElementCollection* vfec;
+
+    //finite element space for pressure
+    mfem::ParFiniteElementSpace* pfes;
+    //finite element collection for pressure
+    mfem::FiniteElementCollection* pfec;
+
+
+    mfem::ConstantCoefficient zeroc;
+    mfem::ConstantCoefficient onec;
+    mfem::Coefficient* mu;
+    mfem::Coefficient* brink; //Brinkman penalization
+
+
+    mfem::ParBilinearForm* af;
+    mfem::ParMixedBilinearForm* bf;
+    std::unique_ptr<mfem::ParBilinearForm> mf;
+
+    mfem::ParLinearForm* lf; //forces
+    mfem::ParLinearForm* pf; //pressure rhs
+
+    // boundary conditions for x,y, and z directions
+    std::map<int, mfem::ConstantCoefficient> bcx;
+    std::map<int, mfem::ConstantCoefficient> bcy;
+    std::map<int, mfem::ConstantCoefficient> bcz;
+
+    // holds BC in coefficient form
+    std::map<int, mfem::Coefficient*> bccx;
+    std::map<int, mfem::Coefficient*> bccy;
+    std::map<int, mfem::Coefficient*> bccz;
+
+    // holds the velocity constrained DOFs
+    mfem::Array<int> ess_tdofv;
+
+    // holds the pressure constrained DOFs
+    mfem::Array<int> ess_tdofp;
+
+    std::unique_ptr<mfem::HypreParMatrix> A;
+    std::unique_ptr<mfem::HypreParMatrix> Ae;
+    std::unique_ptr<mfem::HypreParMatrix> B;
+    std::unique_ptr<mfem::HypreParMatrix> Be;
+    std::unique_ptr<mfem::HypreParMatrix> D;
+    std::unique_ptr<mfem::HypreParMatrix> De;
+    std::unique_ptr<mfem::HypreParMatrix> M;
+    std::unique_ptr<mfem::HypreParMatrix> Me;
+
+    std::unique_ptr<InverseCoeff> ivisc;
+
+    mfem::BlockOperator* bop; //tangent operator
+    mfem::Solver* pop; //preconditioner
+
+    mfem::IterativeSolver *ls;
+
+    void SetEssTDofs(mfem::Vector& bsol, mfem::Array<int>& ess_dofs);
+
+
+    void AssemblePrec1();
+    void AssemblePrec2();
+    void AssemblePrec3();
+    void AssemblePrec4();
+    void AssemblePrec5();
+    void AssemblePrec6();
+    void AssemblePrec7();
+    std::unique_ptr<mfem::HypreBoomerAMG> prec1;
+    std::unique_ptr<mfem::HypreBoomerAMG> prec2;
+
+
+    class LSC:public mfem::Operator
+    {
+    public:
+        LSC(mfem::HypreParMatrix& B_,
+            mfem::HypreParMatrix& A_,
+            mfem::HypreParMatrix& D_,
+            mfem::HypreParMatrix& M_)
+        {
+            B=&B_;
+            A=&A_;
+            D=&D_;
+            M=&M_;
+
+            prec=std::unique_ptr<mfem::HypreBoomerAMG>(new mfem::HypreBoomerAMG());
+            prec->SetOperator(*M);
+            prec->SetPrintLevel(1);
+
+            p.SetSize(M->GetNumCols());
+            u.SetSize(A->GetNumCols());
+            v.SetSize(A->GetNumCols());
+
+            mfem::Operator::width=M->GetNumCols();
+            mfem::Operator::height=M->GetNumCols();
+
+            ls=std::unique_ptr<mfem::IterativeSolver>(new mfem::CGSolver(B->GetComm()));
+            ls->SetOperator(*M);
+            ls->SetPreconditioner(*prec);
+            ls->SetAbsTol(1e-8);
+            ls->SetRelTol(1e-8);
+            ls->SetMaxIter(100);
+            ls->SetPrintLevel(0);
+        }
+
+        virtual
+        ~LSC(){}
+
+        /// Operator application
+        void Mult (const mfem::Vector & x, mfem::Vector & y) const override
+        {
+            prec->Mult(x,p);
+            D->Mult(p,u);
+            A->Mult(u,v);
+            B->Mult(v,p);
+            prec->Mult(p,y);
+        }
+
+        /// Action of the transpose operator
+        void MultTranspose (const mfem::Vector & x, mfem::Vector & y) const override
+        {
+            Mult(x,y);
+        }
+
+    private:
+        mfem::HypreParMatrix *B;
+        mfem::HypreParMatrix *D;
+        mfem::HypreParMatrix *A;
+        mfem::HypreParMatrix *M;
+
+        std::unique_ptr<mfem::HypreBoomerAMG> prec;
+        std::unique_ptr<mfem::IterativeSolver> ls;
+
+        mutable mfem::Vector p;
+        mutable mfem::Vector u;
+        mutable mfem::Vector v;
+    };
+
+
+};
+
+
+
+#endif // STOKESSOLVER_HPP

miniapps/mtop/test_eig.cpp

@@ -0,0 +1,164 @@
+#include "mfem.hpp"
+#include "rand_eigensolver.hpp"
+#include "mtop_solvers.hpp"
+
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI and HYPRE.
+   Mpi::Init();
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../../data/star.mesh";
+   int order = 1;
+   bool static_cond = false;
+   bool pa = false;
+   bool fa = false;
+   const char *device_config = "cpu";
+   bool visualization = true;
+   bool algebraic_ceed = false;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&fa, "-fa", "--full-assembly", "-no-fa",
+                  "--no-full-assembly", "Enable Full Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+#ifdef MFEM_USE_CEED
+   args.AddOption(&algebraic_ceed, "-a", "--algebraic",
+                  "-no-a", "--no-algebraic",
+                  "Use algebraic Ceed solver");
+#endif
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 3. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   if (myid == 0) { device.Print(); }
+
+   // 4. Read the (serial) mesh from the given mesh file on all processors.  We
+   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   //    and volume meshes with the same code.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   // 5. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement. We choose
+   //    'ref_levels' to be the largest number that gives a final mesh with no
+   //    more than 10,000 elements.
+   {
+      int ref_levels =
+         (int)floor(log(10000./mesh.GetNE())/log(2.)/dim);
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh.UniformRefinement();
+      }
+   }
+
+   // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+   {
+      int par_ref_levels = 1;
+      for (int l = 0; l < par_ref_levels; l++)
+      {
+         pmesh.UniformRefinement();
+      }
+   }
+   if(Mpi::WorldRank()==0){
+       std::cout<<pmesh.GetNE()<<std::endl;
+   }
+
+   LElasticOperator* le=new LElasticOperator(&pmesh,1);
+
+   mfem::ConstantCoefficient E(1.0);
+   mfem::ConstantCoefficient nu(0.2);
+   mfem::ConstantCoefficient rho(1.0);
+
+
+   le->AddDispBC(1,4,0.0);
+   le->SetMaterial(E,nu);
+   le->SetVolForce(1.0,1.0);
+   le->Assemble();
+   le->FSolve();
+
+   mfem::ParGridFunction disp;
+   le->GetSol(disp);
+
+   {
+       ParaViewDataCollection paraview_dc("eigp", &pmesh);
+       paraview_dc.SetPrefixPath("ParaView");
+       paraview_dc.SetLevelsOfDetail(order);
+       paraview_dc.SetDataFormat(VTKFormat::BINARY);
+       paraview_dc.SetHighOrderOutput(true);
+       paraview_dc.SetCycle(0);
+       paraview_dc.SetTime(0.0);
+       paraview_dc.RegisterField("disp",&disp);
+       paraview_dc.Save();
+   }
+
+   delete le;
+
+   FRElasticSolver* fr=new FRElasticSolver(&pmesh,2,0.0);
+   fr->AddDispBC(1,4,0.0);
+   fr->SetMaterial(E,nu);
+   fr->SetDensity(rho);
+   fr->Assemble();
+   fr->AssembleSVD();
+
+   const RandomizedSubspaceIteration* ss=fr->GetEigSolver();
+
+   ParaViewDataCollection paraview_dc("eigt", &pmesh);
+   paraview_dc.SetPrefixPath("ParaView");
+   paraview_dc.SetLevelsOfDetail(order);
+   paraview_dc.SetDataFormat(VTKFormat::BINARY);
+   paraview_dc.SetHighOrderOutput(true);
+   paraview_dc.RegisterField("disp",&disp);
+
+   auto eigv=ss->GetModes();
+   fr->GetSol(disp);
+
+   for(int i=0;i<ss->GetNumModes();i++){
+      disp.SetFromTrueDofs(eigv[i]);
+      paraview_dc.SetCycle(i);
+      paraview_dc.SetTime(double(i));
+      paraview_dc.Save();
+   }
+   delete fr;
+
+   Mpi::Finalize();
+   return 0;
+}

miniapps/mtop/test_evec.cpp

@@ -0,0 +1,465 @@
+#include "mfem.hpp"
+#include "mtop_solvers.hpp"
+
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+
+class VectorForceCoeff:public VectorCoefficient
+{
+public:
+    VectorForceCoeff(Vector& A_, Vector& cc, real_t r_)
+        : VectorCoefficient(A_.Size()), a(A_), c(cc), r(r_)
+    {
+    }
+
+    virtual void Eval(Vector &V, ElementTransformation &T,
+                      const IntegrationPoint &ip)
+    {
+        V.SetSize(a.Size());
+        Vector transip(a.Size()); transip=0.0;
+        T.Transform(ip, transip);
+
+
+        transip.Add(-1.0,c);
+        real_t dist=transip.Norml2();
+
+        if(dist<r){
+            V.Set(1.0,a);
+        }else{
+            V=0.0;
+        }
+
+    }
+
+private:
+    Vector a;
+    Vector c;
+    real_t r;
+};
+
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI and HYPRE.
+   Mpi::Init();
+   int num_procs = Mpi::WorldSize();
+   int myrank = Mpi::WorldRank();
+   Hypre::Init();
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../../data/star.mesh";
+   int order = 3;
+   bool static_cond = false;
+   bool pa = false;
+   bool fa = false;
+   const char *device_config = "cpu";
+   bool visualization = true;
+   bool algebraic_ceed = false;
+   real_t freq=0.0;
+   real_t adamp=0.02;
+   real_t bdamp=0.00;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&fa, "-fa", "--full-assembly", "-no-fa",
+                  "--no-full-assembly", "Enable Full Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&freq, "-fr", "--fr",
+                  "Excitation frequency.");
+   args.AddOption(&adamp, "-a", "--alpha",
+                  "alpha*M+beta*K");
+   args.AddOption(&bdamp, "-b", "--beta",
+                  "alpha*M+beta*K");
+#ifdef MFEM_USE_CEED
+   args.AddOption(&algebraic_ceed, "-a", "--algebraic",
+                  "-no-a", "--no-algebraic",
+                  "Use algebraic Ceed solver");
+#endif
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myrank == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myrank == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 3. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   if (myrank == 0) { device.Print(); }
+
+   // 4. Read the (serial) mesh from the given mesh file on all processors.  We
+   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   //    and volume meshes with the same code.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   // 5. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement. We choose
+   //    'ref_levels' to be the largest number that gives a final mesh with no
+   //    more than 10,000 elements.
+   {
+      int ref_levels =
+         (int)floor(log(1000./mesh.GetNE())/log(2.)/dim);
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh.UniformRefinement();
+      }
+   }
+
+   // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+   {
+      int par_ref_levels = 0;
+      for (int l = 0; l < par_ref_levels; l++)
+      {
+         pmesh.UniformRefinement();
+      }
+   }
+   if(Mpi::WorldRank()==0){
+       std::cout<<pmesh.GetNE()<<std::endl;
+   }
+
+   H1_FECollection* vfec=new H1_FECollection(order,dim);
+   //ParFiniteElementSpace* vfes=new ParFiniteElementSpace(&pmesh,vfec,dim, mfem::Ordering::byVDIM);
+   ParFiniteElementSpace* vfes=new ParFiniteElementSpace(&pmesh,vfec,dim, mfem::Ordering::byNODES);
+
+   //extract the BC dofs
+   Array<int> ess_dofs; //true dofs
+   {
+       Array<int> ess_bdr(pmesh.bdr_attributes.Max());
+       ess_bdr = 0;
+       ess_bdr[0] = 1;
+       ess_bdr[1] = 1;
+       ess_bdr[2] = 1;
+       ess_bdr[3] = 1;
+       ess_bdr[4] = 0;
+       vfes->GetEssentialTrueDofs(ess_bdr, ess_dofs);
+   }
+
+
+   std::unique_ptr<ParBilinearForm> mf; //mass bilinear form
+   std::unique_ptr<ParBilinearForm> cf; //damping bilinear form
+   std::unique_ptr<ParBilinearForm> kf; //stiffness bilinear form
+   std::unique_ptr<ParBilinearForm> wf;
+
+   ConstantCoefficient rho(1.0);
+   ConstantCoefficient damp(adamp);
+   ConstantCoefficient E(1.0);
+   ConstantCoefficient nu(0.2);
+
+   IsoElasticyLambdaCoeff lambda(&E,&nu);
+   IsoElasticySchearCoeff mu(&E,&nu);
+
+   //stiffness operator
+   kf.reset(new ParBilinearForm(vfes));
+   //mass operator
+   mf.reset(new ParBilinearForm(vfes));
+   //damping operator
+   cf.reset(new ParBilinearForm(vfes));
+   //Helmholtz operator
+   wf.reset(new ParBilinearForm(vfes));
+
+   if(pa){
+       kf->SetAssemblyLevel(mfem::AssemblyLevel::PARTIAL);
+       mf->SetAssemblyLevel(mfem::AssemblyLevel::PARTIAL);
+       cf->SetAssemblyLevel(mfem::AssemblyLevel::PARTIAL);
+       wf->SetAssemblyLevel(mfem::AssemblyLevel::PARTIAL);
+   }else{
+       kf->SetAssemblyLevel(mfem::AssemblyLevel::LEGACY);
+       mf->SetAssemblyLevel(mfem::AssemblyLevel::LEGACY);
+       cf->SetAssemblyLevel(mfem::AssemblyLevel::LEGACY);
+       wf->SetAssemblyLevel(mfem::AssemblyLevel::LEGACY);
+   }
+
+   std::cout<<"Add Integrators"<<std::endl;
+
+   kf->AddDomainIntegrator(new ElasticityIntegrator(lambda,mu));
+   mf->AddDomainIntegrator(new VectorMassIntegrator(rho));
+   cf->AddDomainIntegrator(new VectorMassIntegrator(damp));
+
+   ConstantCoefficient bE(bdamp);
+   ConstantCoefficient bnu(0.2);
+
+   IsoElasticyLambdaCoeff blambda(&bE,&bnu);
+   IsoElasticySchearCoeff bmu(&bE,&bnu);
+   cf->AddDomainIntegrator(new ElasticityIntegrator(blambda,bmu));
+   //Helmholtz operator
+   wf->AddDomainIntegrator(new ElasticityIntegrator(lambda,mu));
+   ProductCoefficient pc(-freq*freq,rho);
+   wf->AddDomainIntegrator(new VectorMassIntegrator(pc));
+
+
+   std::cout<<"Start Assemble"<<std::endl;
+   kf->Assemble();
+   kf->Finalize();
+   std::cout<<"K ready"<<std::endl;
+   mf->Assemble();
+   mf->Finalize();
+   std::cout<<"M ready"<<std::endl;
+   cf->Assemble();
+   cf->Finalize();
+   std::cout<<"C ready"<<std::endl;
+   wf->Assemble();
+   wf->Finalize();
+   std::cout<<"W ready"<<std::endl;
+
+   std::cout<<"Cons"<<std::endl;
+
+   //(-M*\omega^2+i*\omega*C+K)u=f
+
+   unique_ptr<HypreParMatrix> kmat; kmat.reset(kf->ParallelAssemble());
+   unique_ptr<HypreParMatrix> cmat; cmat.reset(cf->ParallelAssemble());
+   unique_ptr<HypreParMatrix> mmat; mmat.reset(mf->ParallelAssemble());
+   unique_ptr<HypreParMatrix> wmat; wmat.reset(wf->ParallelAssemble());
+
+   kmat->EliminateBC(ess_dofs,Operator::DiagonalPolicy::DIAG_ONE);
+   mmat->EliminateBC(ess_dofs,Operator::DiagonalPolicy::DIAG_ZERO);
+   cmat->EliminateBC(ess_dofs,Operator::DiagonalPolicy::DIAG_ONE);
+   wmat->EliminateBC(ess_dofs,Operator::DiagonalPolicy::DIAG_ONE);
+
+   //kmat->EliminateRowsCols(ess_dofs); kmat->EliminateZeroRows();
+   //mmat->EliminateRowsCols(ess_dofs);
+   //cmat->EliminateRowsCols(ess_dofs);
+   //wmat->EliminateRowsCols(ess_dofs); wmat->EliminateZeroRows();
+
+
+   const Operator *P = kf->GetProlongation();
+   ConstrainedOperator ckf(new RAPOperator(*P,*kf,*P),ess_dofs,true, Operator::DiagonalPolicy::DIAG_ONE);
+   ConstrainedOperator cmf(new RAPOperator(*P,*mf,*P),ess_dofs,true, Operator::DiagonalPolicy::DIAG_ZERO);
+   ConstrainedOperator ccf(new RAPOperator(*P,*cf,*P),ess_dofs,true, Operator::DiagonalPolicy::DIAG_ONE);
+   ConstrainedOperator cwf(new RAPOperator(*P,*wf,*P),ess_dofs,true, Operator::DiagonalPolicy::DIAG_ONE);
+
+   OperatorHandle hkf;
+   OperatorHandle hmf;
+   OperatorHandle hcf;
+   OperatorHandle hwf;
+   kf->FormSystemMatrix(ess_dofs, hkf);
+   mf->FormSystemMatrix(ess_dofs, hmf);
+   cf->FormSystemMatrix(ess_dofs, hcf);
+   wf->FormSystemMatrix(ess_dofs, hwf);
+
+
+   unique_ptr<HypreParMatrix> bm;
+   {
+       Array2D<const HypreParMatrix*> am(2,2);
+       am(0,0)=wmat.get();
+       am(0,1)=cmat.get();
+       am(1,0)=cmat.get();
+       am(1,1)=wmat.get();
+
+       Array2D<real_t> cm(2,2);
+       cm(0,0)=1.0; cm(0,1)=-1;
+       cm(1,0)=1.0; cm(1,1)=1;
+
+       bm.reset(HypreParMatrixFromBlocks(am,&cm));
+    }
+
+   std::cout<<"Forcing!"<<std::endl;
+
+   //force coefficient
+   Vector A(dim); A=0.0; A[1]=1.0; //force amplitude
+   Vector cc(dim); cc[0]=3.9; cc[1]=0.5; if(dim==3){cc[2]=0.0;} //force possition
+   VectorForceCoeff fc(A,cc,0.1);
+
+   ParLinearForm lf(vfes);
+   lf.AddDomainIntegrator(new VectorDomainLFIntegrator(fc));
+   lf.Assemble();
+
+   std::cout<<"LF is assmebled!"<<std::endl;
+
+   Array<int> block_true_offsets;
+   block_true_offsets.SetSize(3);
+   block_true_offsets[0] = 0;
+   block_true_offsets[1] = vfes->GetTrueVSize();
+   block_true_offsets[2] = vfes->GetTrueVSize();
+   block_true_offsets.PartialSum();
+
+   // solution
+   BlockVector x; x.Update(block_true_offsets); x=0.0;
+   // RHS
+   BlockVector f; f.Update(block_true_offsets);
+   std::cout<<pmesh.GetMyRank()<<" f0.size="<<f.GetBlock(0).Size()<<std::endl;
+   std::cout<<pmesh.GetMyRank()<<" f1.size="<<f.GetBlock(1).Size()<<std::endl;
+   std::cout<<pmesh.GetMyRank()<<" cwf.size="<<cwf.Width()<<std::endl;
+   std::cout<<pmesh.GetMyRank()<<" vfes.size="<<vfes->GetTrueVSize()<<std::endl;
+   lf.ParallelAssemble(f.GetBlock(0));
+   lf.ParallelAssemble(f.GetBlock(1)); f.GetBlock(1)=0.0;
+
+   //mass proportional force
+   x.GetBlock(0)=1.0;
+   mmat->Mult(x.GetBlock(0),f.GetBlock(0));
+   x=0.0;
+
+
+   std::cout<<"RHS is ready!"<<std::endl;
+
+   //cwf.EliminateRHS(x.GetBlock(0), f.GetBlock(0));
+   //cwf.EliminateRHS(x.GetBlock(1), f.GetBlock(1));
+
+   BlockOperator bop(block_true_offsets);
+   bop.SetBlock(0,0,wmat.get(),1.0);
+   bop.SetBlock(0,1,cmat.get(),-1.0);
+   bop.SetBlock(1,0,cmat.get(),1.0);
+   bop.SetBlock(1,1,wmat.get(),1.0);
+
+
+   //get the first n eigenmodes
+   MUMPSSolver* mumps=new MUMPSSolver(pmesh.GetComm());
+   mumps->SetPrintLevel(1);
+   mumps->SetMatrixSymType(MUMPSSolver::MatType::SYMMETRIC_POSITIVE_DEFINITE);
+   mumps->SetOperator(*kmat);
+
+   int num_evec=20;
+   AdaptiveRandomizedGenEig ae(pmesh.GetComm());
+   ae.SetOperators(*mmat,*kmat,*mumps);
+   ae.SetNumModes(num_evec+10);
+   ae.SetNumIter(4);
+   ae.SolveNA();
+
+   delete mumps;
+
+   std::cout<<"Allocate EVECP"<<std::endl;
+
+   EVECPrec* prec1=new EVECPrec(pmesh.GetComm());
+   prec1->SetOperators(kmat.get(),mmat.get(),cmat.get(), ae.GetModes(), num_evec, 1.0, freq*freq, 1.0);
+
+
+   MSP1Prec* prec2=new MSP1Prec(pmesh.GetComm());
+   prec2->SetOperators(kmat.get(),mmat.get(),cmat.get(),1.0,freq*freq,1.0);
+   prec2->SetAbsTol(1e-12);
+   prec2->SetRelTol(1e-6);
+   prec2->SetMaxIter(1);
+
+   ForwardOp* prec3=new ForwardOp(pmesh.GetComm());
+   prec3->SetOperators(kmat.get(),mmat.get(),cmat.get(),1.0,freq*freq,1.0);
+
+   ProductSolver* prec4=new ProductSolver(&bop,prec1,prec2,false,false,false);
+
+
+
+   Solver* prec=prec4;
+
+   prec->Mult(f,x);
+   {
+       Vector res(f.Size()); res=0.0;
+       bop.Mult(x,res);
+       res.Add(-1,f);
+       real_t rr=mfem::InnerProduct(pmesh.GetComm(), res,res);
+       real_t fr=mfem::InnerProduct(pmesh.GetComm(), f,f);
+       if(pmesh.GetMyRank()==0){
+           std::cout<<"Residual="<<sqrt(rr)<<" |f|="<<sqrt(fr)<<std::endl;
+           std::cout<<"fr="<<freq<<" fr^2="<<freq*freq<<std::endl;
+       }
+   }
+
+
+   Vector xstat(f.GetBlock(0));xstat=0.0;
+
+
+
+   /*
+   delete prec;
+   delete vfes;
+   delete vfec;
+   Mpi::Finalize();
+   return 0;
+   */
+
+
+   //set the linear solver
+   FGMRESSolver* ls=new FGMRESSolver(pmesh.GetComm());
+   //ls->SetOperator(bop);
+   ls->SetOperator(*bm);
+   ls->SetAbsTol(1e-12);
+   ls->SetRelTol(1e-12);
+   ls->SetMaxIter(1000);
+   ls->SetPrintLevel(1);
+   ls->SetKDim(1000);
+
+   ls->SetOperator(bop);
+   ls->SetPreconditioner(*prec);
+   ls->Mult(f,x);
+   delete ls;
+
+   delete prec1;
+   delete prec2;
+   delete prec3;
+   delete prec4;
+
+   //check the solution
+   {
+       Vector xm(x); xm=0.0;
+       MUMPSSolver mumps(bm->GetComm());
+       mumps.SetPrintLevel(2);
+       mumps.SetMatrixSymType(MUMPSSolver::MatType::UNSYMMETRIC);
+       mumps.SetOperator(*bm);
+       mumps.Mult(f, xm);
+
+       xm.Add(-1,x);
+       real_t rr=mfem::InnerProduct(pmesh.GetComm(), xm,xm);
+       if(pmesh.GetMyRank()==0){
+           std::cout<<"|xm-x|="<<sqrt(rr)<<std::endl;
+           std::cout<<"fr="<<freq<<"       fr^2="<<freq*freq<<std::endl;
+       }
+   }
+
+
+   //check the solutions
+   ParGridFunction rx(vfes); rx.SetFromTrueDofs(x.GetBlock(0));
+   ParGridFunction ix(vfes); ix.SetFromTrueDofs(x.GetBlock(1)); ix*=-1.0;
+   ParGridFunction xs(vfes); xs.SetFromTrueDofs(xstat);
+   {
+        ParaViewDataCollection paraview_dc("stokes", &pmesh);
+        paraview_dc.SetPrefixPath("ParaView");
+        paraview_dc.SetLevelsOfDetail(order);
+        paraview_dc.SetDataFormat(VTKFormat::BINARY);
+        paraview_dc.SetHighOrderOutput(true);
+        paraview_dc.SetCycle(0);
+        paraview_dc.SetTime(0.0);
+        paraview_dc.RegisterField("re",&rx);
+        paraview_dc.RegisterField("im",&ix);
+        paraview_dc.RegisterField("xs",&xs);
+        paraview_dc.Save();
+   }
+
+   BlockVector r(f);
+
+   bop.Mult(x,r);
+   r.Add(-1,f);
+   real_t rr=mfem::InnerProduct(pmesh.GetComm(), r,r);
+   if(pmesh.GetMyRank()==0){
+       std::cout<<"Residual="<<sqrt(rr)<<std::endl;
+       std::cout<<"fr="<<freq<<"       fr^2="<<freq*freq<<std::endl;
+   }
+
+   delete vfes;
+   delete vfec;
+   Mpi::Finalize();
+   return 0;
+}

miniapps/mtop/test_msp1.cpp

@@ -0,0 +1,439 @@
+#include "mfem.hpp"
+#include "mtop_solvers.hpp"
+
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+
+class VectorForceCoeff:public VectorCoefficient
+{
+public:
+    VectorForceCoeff(Vector& A_, Vector& cc, real_t r_)
+        : VectorCoefficient(A_.Size()), a(A_), c(cc), r(r_)
+    {
+    }
+
+    virtual void Eval(Vector &V, ElementTransformation &T,
+                      const IntegrationPoint &ip)
+    {
+        V.SetSize(a.Size());
+        Vector transip(a.Size()); transip=0.0;
+        T.Transform(ip, transip);
+
+
+        transip.Add(-1.0,c);
+        real_t dist=transip.Norml2();
+
+        if(dist<r){
+            V.Set(1.0,a);
+        }else{
+            V=0.0;
+        }
+
+    }
+
+private:
+    Vector a;
+    Vector c;
+    real_t r;
+};
+
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI and HYPRE.
+   Mpi::Init();
+   int num_procs = Mpi::WorldSize();
+   int myrank = Mpi::WorldRank();
+   Hypre::Init();
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../../data/star.mesh";
+   int order = 3;
+   bool static_cond = false;
+   bool pa = false;
+   bool fa = false;
+   const char *device_config = "cpu";
+   bool visualization = true;
+   bool algebraic_ceed = false;
+   real_t freq=0.0;
+   real_t adamp=0.02;
+   real_t bdamp=0.00;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&fa, "-fa", "--full-assembly", "-no-fa",
+                  "--no-full-assembly", "Enable Full Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&freq, "-fr", "--fr",
+                  "Excitation frequency.");
+   args.AddOption(&adamp, "-a", "--alpha",
+                  "alpha*M+beta*K");
+   args.AddOption(&bdamp, "-b", "--beta",
+                  "alpha*M+beta*K");
+#ifdef MFEM_USE_CEED
+   args.AddOption(&algebraic_ceed, "-a", "--algebraic",
+                  "-no-a", "--no-algebraic",
+                  "Use algebraic Ceed solver");
+#endif
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myrank == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myrank == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 3. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   if (myrank == 0) { device.Print(); }
+
+   // 4. Read the (serial) mesh from the given mesh file on all processors.  We
+   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   //    and volume meshes with the same code.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   // 5. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement. We choose
+   //    'ref_levels' to be the largest number that gives a final mesh with no
+   //    more than 10,000 elements.
+   {
+      int ref_levels =
+         (int)floor(log(1000./mesh.GetNE())/log(2.)/dim);
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh.UniformRefinement();
+      }
+   }
+
+   // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+   {
+      int par_ref_levels = 0;
+      for (int l = 0; l < par_ref_levels; l++)
+      {
+         pmesh.UniformRefinement();
+      }
+   }
+   if(Mpi::WorldRank()==0){
+       std::cout<<pmesh.GetNE()<<std::endl;
+   }
+
+   H1_FECollection* vfec=new H1_FECollection(order,dim);
+   //ParFiniteElementSpace* vfes=new ParFiniteElementSpace(&pmesh,vfec,dim, mfem::Ordering::byVDIM);
+   ParFiniteElementSpace* vfes=new ParFiniteElementSpace(&pmesh,vfec,dim, mfem::Ordering::byNODES);
+
+   //extract the BC dofs
+   Array<int> ess_dofs; //true dofs
+   {
+       Array<int> ess_bdr(pmesh.bdr_attributes.Max());
+       ess_bdr = 0;
+       ess_bdr[0] = 1;
+       ess_bdr[1] = 1;
+       ess_bdr[2] = 1;
+       ess_bdr[3] = 1;
+       ess_bdr[4] = 0;
+       vfes->GetEssentialTrueDofs(ess_bdr, ess_dofs);
+   }
+
+
+   std::unique_ptr<ParBilinearForm> mf; //mass bilinear form
+   std::unique_ptr<ParBilinearForm> cf; //damping bilinear form
+   std::unique_ptr<ParBilinearForm> kf; //stiffness bilinear form
+   std::unique_ptr<ParBilinearForm> wf;
+
+   ConstantCoefficient rho(1.0);
+   ConstantCoefficient damp(adamp);
+   ConstantCoefficient E(1.0);
+   ConstantCoefficient nu(0.2);
+
+   IsoElasticyLambdaCoeff lambda(&E,&nu);
+   IsoElasticySchearCoeff mu(&E,&nu);
+
+   //stiffness operator
+   kf.reset(new ParBilinearForm(vfes));
+   //mass operator
+   mf.reset(new ParBilinearForm(vfes));
+   //damping operator
+   cf.reset(new ParBilinearForm(vfes));
+   //Helmholtz operator
+   wf.reset(new ParBilinearForm(vfes));
+
+   if(pa){
+       kf->SetAssemblyLevel(mfem::AssemblyLevel::PARTIAL);
+       mf->SetAssemblyLevel(mfem::AssemblyLevel::PARTIAL);
+       cf->SetAssemblyLevel(mfem::AssemblyLevel::PARTIAL);
+       wf->SetAssemblyLevel(mfem::AssemblyLevel::PARTIAL);
+   }else{
+       kf->SetAssemblyLevel(mfem::AssemblyLevel::LEGACY);
+       mf->SetAssemblyLevel(mfem::AssemblyLevel::LEGACY);
+       cf->SetAssemblyLevel(mfem::AssemblyLevel::LEGACY);
+       wf->SetAssemblyLevel(mfem::AssemblyLevel::LEGACY);
+   }
+
+   std::cout<<"Add Integrators"<<std::endl;
+
+   kf->AddDomainIntegrator(new ElasticityIntegrator(lambda,mu));
+   mf->AddDomainIntegrator(new VectorMassIntegrator(rho));
+   cf->AddDomainIntegrator(new VectorMassIntegrator(damp));
+
+   ConstantCoefficient bE(bdamp);
+   ConstantCoefficient bnu(0.2);
+
+   IsoElasticyLambdaCoeff blambda(&bE,&bnu);
+   IsoElasticySchearCoeff bmu(&bE,&bnu);
+   cf->AddDomainIntegrator(new ElasticityIntegrator(blambda,bmu));
+   //Helmholtz operator
+   wf->AddDomainIntegrator(new ElasticityIntegrator(lambda,mu));
+   ProductCoefficient pc(-freq*freq,rho);
+   //ConstantCoefficient pc(-freq*freq*1.0);
+   wf->AddDomainIntegrator(new VectorMassIntegrator(pc));
+
+
+   std::cout<<"Start Assemble"<<std::endl;
+   kf->Assemble();
+   kf->Finalize();
+   std::cout<<"K ready"<<std::endl;
+   mf->Assemble();
+   mf->Finalize();
+   std::cout<<"M ready"<<std::endl;
+   cf->Assemble();
+   cf->Finalize();
+   std::cout<<"C ready"<<std::endl;
+   wf->Assemble();
+   wf->Finalize();
+   std::cout<<"W ready"<<std::endl;
+
+   std::cout<<"Cons"<<std::endl;
+
+   //(-M*\omega^2+i*\omega*C+K)u=f
+
+   unique_ptr<HypreParMatrix> kmat; kmat.reset(kf->ParallelAssemble());
+   unique_ptr<HypreParMatrix> cmat; cmat.reset(cf->ParallelAssemble());
+   unique_ptr<HypreParMatrix> mmat; mmat.reset(mf->ParallelAssemble());
+   unique_ptr<HypreParMatrix> wmat; wmat.reset(wf->ParallelAssemble());
+
+   kmat->EliminateBC(ess_dofs,Operator::DiagonalPolicy::DIAG_ONE);
+   mmat->EliminateBC(ess_dofs,Operator::DiagonalPolicy::DIAG_ZERO);
+   cmat->EliminateBC(ess_dofs,Operator::DiagonalPolicy::DIAG_ONE);
+   wmat->EliminateBC(ess_dofs,Operator::DiagonalPolicy::DIAG_ONE);
+
+   //kmat->EliminateRowsCols(ess_dofs); kmat->EliminateZeroRows();
+   //mmat->EliminateRowsCols(ess_dofs);
+   //cmat->EliminateRowsCols(ess_dofs);
+   //wmat->EliminateRowsCols(ess_dofs); wmat->EliminateZeroRows();
+
+
+   const Operator *P = kf->GetProlongation();
+   ConstrainedOperator ckf(new RAPOperator(*P,*kf,*P),ess_dofs,true, Operator::DiagonalPolicy::DIAG_ONE);
+   ConstrainedOperator cmf(new RAPOperator(*P,*mf,*P),ess_dofs,true, Operator::DiagonalPolicy::DIAG_ZERO);
+   ConstrainedOperator ccf(new RAPOperator(*P,*cf,*P),ess_dofs,true, Operator::DiagonalPolicy::DIAG_ONE);
+   ConstrainedOperator cwf(new RAPOperator(*P,*wf,*P),ess_dofs,true, Operator::DiagonalPolicy::DIAG_ONE);
+
+   OperatorHandle hkf;
+   OperatorHandle hmf;
+   OperatorHandle hcf;
+   OperatorHandle hwf;
+   kf->FormSystemMatrix(ess_dofs, hkf);
+   mf->FormSystemMatrix(ess_dofs, hmf);
+   cf->FormSystemMatrix(ess_dofs, hcf);
+   wf->FormSystemMatrix(ess_dofs, hwf);
+
+
+   unique_ptr<HypreParMatrix> bm;
+   {
+       Array2D<const HypreParMatrix*> am(2,2);
+       am(0,0)=wmat.get();
+       am(0,1)=cmat.get();
+       am(1,0)=cmat.get();
+       am(1,1)=wmat.get();
+
+       Array2D<real_t> cm(2,2);
+       cm(0,0)=1.0; cm(0,1)=-1;
+       cm(1,0)=1.0; cm(1,1)=1;
+
+       bm.reset(HypreParMatrixFromBlocks(am,&cm));
+    }
+
+   std::cout<<"Forcing!"<<std::endl;
+
+   //force coefficient
+   Vector A(dim); A=0.0; A[1]=1.0; //force amplitude
+   Vector cc(dim); cc[0]=3.9; cc[1]=0.5; if(dim==3){cc[2]=0.0;} //force possition
+   VectorForceCoeff fc(A,cc,0.1);
+
+   ParLinearForm lf(vfes);
+   lf.AddDomainIntegrator(new VectorDomainLFIntegrator(fc));
+   lf.Assemble();
+
+   std::cout<<"LF is assmebled!"<<std::endl;
+
+   Array<int> block_true_offsets;
+   block_true_offsets.SetSize(3);
+   block_true_offsets[0] = 0;
+   block_true_offsets[1] = vfes->GetTrueVSize();
+   block_true_offsets[2] = vfes->GetTrueVSize();
+   block_true_offsets.PartialSum();
+
+   // solution
+   BlockVector x; x.Update(block_true_offsets); x=0.0;
+   // RHS
+   BlockVector f; f.Update(block_true_offsets);
+   std::cout<<pmesh.GetMyRank()<<" f0.size="<<f.GetBlock(0).Size()<<std::endl;
+   std::cout<<pmesh.GetMyRank()<<" f1.size="<<f.GetBlock(1).Size()<<std::endl;
+   std::cout<<pmesh.GetMyRank()<<" cwf.size="<<cwf.Width()<<std::endl;
+   std::cout<<pmesh.GetMyRank()<<" vfes.size="<<vfes->GetTrueVSize()<<std::endl;
+   lf.ParallelAssemble(f.GetBlock(0));
+   lf.ParallelAssemble(f.GetBlock(1)); f.GetBlock(1)=0.0;
+
+   std::cout<<"RHS is ready!"<<std::endl;
+
+   //cwf.EliminateRHS(x.GetBlock(0), f.GetBlock(0));
+   //cwf.EliminateRHS(x.GetBlock(1), f.GetBlock(1));
+
+   BlockOperator bop(block_true_offsets);
+   bop.SetBlock(0,0,wmat.get(),1.0);
+   bop.SetBlock(0,1,cmat.get(),-1.0);
+   bop.SetBlock(1,0,cmat.get(),1.0);
+   bop.SetBlock(1,1,wmat.get(),1.0);
+
+
+   SumOperator W(mmat.get(),1.0,mmat.get(),0.0,false,false);
+
+   ParLORDiscretization lork(*kf,ess_dofs);
+   HypreParMatrix& lorkm=lork.GetAssembledMatrix();
+
+   ParFiniteElementSpace& lorfes=lork.GetParFESpace();
+
+   ParLORDiscretization lorc(*cf,ess_dofs);
+   HypreParMatrix& lorcm=lorc.GetAssembledMatrix();
+
+
+   std::cout<<"Allocate MSP1"<<std::endl;
+   /*
+   MSP1Prec* prec=new MSP1Prec(pmesh.GetComm());
+   prec->SetOperators(kmat.get(),mmat.get(),cmat.get(),1.0,freq*freq,1.0);
+   */
+
+   MSP3Prec* prec=new MSP3Prec(pmesh.GetComm());
+   prec->SetOperators(kmat.get(),mmat.get(),cmat.get(),1.0,freq*freq,1.0,10000.0);
+   //prec->SetOperators(&lorkm,&lorcm,1.0,1.0);
+   prec->SetAbsTol(1e-12);
+   prec->SetRelTol(1e-6);
+   prec->SetMaxIter(1);
+
+
+   prec->Mult(f,x);
+   {
+       real_t rr=mfem::InnerProduct(pmesh.GetComm(), x,x);
+       if(pmesh.GetMyRank()==0){
+           std::cout<<"Residual="<<sqrt(rr)<<std::endl;
+           std::cout<<"fr="<<freq<<" fr^2="<<freq*freq<<std::endl;
+       }
+   }
+
+
+   Vector xstat(f.GetBlock(0));xstat=0.0;
+
+   /*
+   delete prec;
+   delete vfes;
+   delete vfec;
+   Mpi::Finalize();
+   return 0;
+    */
+
+   //set the linear solver
+   FGMRESSolver* ls=new FGMRESSolver(pmesh.GetComm());
+   //ls->SetOperator(bop);
+   ls->SetOperator(*bm);
+   ls->SetAbsTol(1e-12);
+   ls->SetRelTol(1e-12);
+   ls->SetMaxIter(1000);
+   ls->SetPrintLevel(1);
+   ls->SetKDim(1000);
+
+   ls->SetOperator(bop);
+   ls->SetPreconditioner(*prec);
+   ls->Mult(f,x);
+   delete ls;
+
+   delete prec;
+
+   //check the solution
+   {
+       Vector xm(x); xm=0.0;
+       MUMPSSolver mumps(bm->GetComm());
+       mumps.SetPrintLevel(2);
+       mumps.SetMatrixSymType(MUMPSSolver::MatType::UNSYMMETRIC);
+       mumps.SetOperator(*bm);
+       mumps.Mult(f, xm);
+
+       xm.Add(-1,x);
+       real_t rr=mfem::InnerProduct(pmesh.GetComm(), xm,xm);
+       if(pmesh.GetMyRank()==0){
+           std::cout<<"|xm-x|="<<sqrt(rr)<<std::endl;
+           std::cout<<"fr="<<freq<<"       fr^2="<<freq*freq<<std::endl;
+       }
+   }
+
+
+   //check the solutions
+   ParGridFunction rx(vfes); rx.SetFromTrueDofs(x.GetBlock(0));
+   ParGridFunction ix(vfes); ix.SetFromTrueDofs(x.GetBlock(1)); ix*=-1.0;
+   ParGridFunction xs(vfes); xs.SetFromTrueDofs(xstat);
+   {
+        ParaViewDataCollection paraview_dc("stokes", &pmesh);
+        paraview_dc.SetPrefixPath("ParaView");
+        paraview_dc.SetLevelsOfDetail(order);
+        paraview_dc.SetDataFormat(VTKFormat::BINARY);
+        paraview_dc.SetHighOrderOutput(true);
+        paraview_dc.SetCycle(0);
+        paraview_dc.SetTime(0.0);
+        paraview_dc.RegisterField("re",&rx);
+        paraview_dc.RegisterField("im",&ix);
+        paraview_dc.RegisterField("xs",&xs);
+        paraview_dc.Save();
+   }
+
+   BlockVector r(f);
+
+   bop.Mult(x,r);
+   r.Add(-1,f);
+   real_t rr=mfem::InnerProduct(pmesh.GetComm(), r,r);
+   if(pmesh.GetMyRank()==0){
+       std::cout<<"Residual="<<sqrt(rr)<<std::endl;
+       std::cout<<"fr="<<freq<<"       fr^2="<<freq*freq<<std::endl;
+   }
+
+   delete vfes;
+   delete vfec;
+   Mpi::Finalize();
+   return 0;
+}

miniapps/mtop/test_presb.cpp

@@ -0,0 +1,638 @@
+#include "mfem.hpp"
+#include "mtop_solvers.hpp"
+
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+
+class VectorForceCoeff:public VectorCoefficient
+{
+public:
+    VectorForceCoeff(Vector& A_, Vector& cc, real_t r_)
+        : VectorCoefficient(A_.Size()), a(A_), c(cc), r(r_)
+    {
+    }
+
+    virtual void Eval(Vector &V, ElementTransformation &T,
+                      const IntegrationPoint &ip)
+    {
+        V.SetSize(a.Size());
+        Vector transip(a.Size()); transip=0.0;
+        T.Transform(ip, transip);
+
+
+        transip.Add(-1.0,c);
+        real_t dist=transip.Norml2();
+
+        if(dist<r){
+            V.Set(1.0,a);
+        }else{
+            V=0.0;
+        }
+
+    }
+
+private:
+    Vector a;
+    Vector c;
+    real_t r;
+};
+
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI and HYPRE.
+   Mpi::Init();
+   int num_procs = Mpi::WorldSize();
+   int myrank = Mpi::WorldRank();
+   Hypre::Init();
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../../data/star.mesh";
+   int order = 3;
+   bool static_cond = false;
+   bool pa = false;
+   bool fa = false;
+   const char *device_config = "cpu";
+   bool visualization = true;
+   bool algebraic_ceed = false;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&fa, "-fa", "--full-assembly", "-no-fa",
+                  "--no-full-assembly", "Enable Full Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+#ifdef MFEM_USE_CEED
+   args.AddOption(&algebraic_ceed, "-a", "--algebraic",
+                  "-no-a", "--no-algebraic",
+                  "Use algebraic Ceed solver");
+#endif
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myrank == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myrank == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 3. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   if (myrank == 0) { device.Print(); }
+
+   // 4. Read the (serial) mesh from the given mesh file on all processors.  We
+   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   //    and volume meshes with the same code.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   // 5. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement. We choose
+   //    'ref_levels' to be the largest number that gives a final mesh with no
+   //    more than 10,000 elements.
+   {
+      int ref_levels =
+         (int)floor(log(1000./mesh.GetNE())/log(2.)/dim);
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh.UniformRefinement();
+      }
+   }
+
+   // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+   {
+      int par_ref_levels = 0;
+      for (int l = 0; l < par_ref_levels; l++)
+      {
+         pmesh.UniformRefinement();
+      }
+   }
+   if(Mpi::WorldRank()==0){
+       std::cout<<pmesh.GetNE()<<std::endl;
+   }
+
+   H1_FECollection* vfec=new H1_FECollection(order,dim);
+   //ParFiniteElementSpace* vfes=new ParFiniteElementSpace(&pmesh,vfec,dim, mfem::Ordering::byVDIM);
+   ParFiniteElementSpace* vfes=new ParFiniteElementSpace(&pmesh,vfec,dim, mfem::Ordering::byNODES);
+
+   //extract the BC dofs
+   Array<int> ess_dofs; //true dofs
+   {
+       Array<int> ess_bdr(pmesh.bdr_attributes.Max());
+       ess_bdr = 0;
+       ess_bdr[0] = 1;
+       ess_bdr[1] = 1;
+       ess_bdr[2] = 1;
+       ess_bdr[3] = 1;
+       ess_bdr[4] = 0;
+       vfes->GetEssentialTrueDofs(ess_bdr, ess_dofs);
+   }
+
+
+   std::unique_ptr<ParBilinearForm> mf; //mass bilinear form
+   std::unique_ptr<ParBilinearForm> cf; //damping bilinear form
+   std::unique_ptr<ParBilinearForm> kf; //stiffness bilinear form
+   std::unique_ptr<ParBilinearForm> wf;
+
+   ConstantCoefficient rho(1.0);
+   ConstantCoefficient damp(0.02);
+   ConstantCoefficient E(1.0);
+   ConstantCoefficient nu(0.2);
+
+   IsoElasticyLambdaCoeff lambda(&E,&nu);
+   IsoElasticySchearCoeff mu(&E,&nu);
+
+   //stiffness operator
+   kf.reset(new ParBilinearForm(vfes));
+   //mass operator
+   mf.reset(new ParBilinearForm(vfes));
+   //damping operator
+   cf.reset(new ParBilinearForm(vfes));
+   //Helmholtz operator
+   wf.reset(new ParBilinearForm(vfes));
+
+   if(pa){
+       kf->SetAssemblyLevel(mfem::AssemblyLevel::PARTIAL);
+       mf->SetAssemblyLevel(mfem::AssemblyLevel::PARTIAL);
+       cf->SetAssemblyLevel(mfem::AssemblyLevel::PARTIAL);
+       wf->SetAssemblyLevel(mfem::AssemblyLevel::PARTIAL);
+   }else{
+       kf->SetAssemblyLevel(mfem::AssemblyLevel::LEGACY);
+       mf->SetAssemblyLevel(mfem::AssemblyLevel::LEGACY);
+       cf->SetAssemblyLevel(mfem::AssemblyLevel::LEGACY);
+       wf->SetAssemblyLevel(mfem::AssemblyLevel::LEGACY);
+   }
+
+   std::cout<<"Add Integrators"<<std::endl;
+
+   kf->AddDomainIntegrator(new ElasticityIntegrator(lambda,mu));
+   mf->AddDomainIntegrator(new VectorMassIntegrator(rho));
+   cf->AddDomainIntegrator(new VectorMassIntegrator(damp));
+   //Helmholtz operator
+   wf->AddDomainIntegrator(new ElasticityIntegrator(lambda,mu));
+   real_t freq=0.5;
+   //ProductCoefficient pc(-freq*freq,rho);
+   ConstantCoefficient pc(-freq*freq*1.0);
+   wf->AddDomainIntegrator(new VectorMassIntegrator(pc));
+
+
+   std::cout<<"Start Assemble"<<std::endl;
+   kf->Assemble();
+   kf->Finalize();
+   std::cout<<"K ready"<<std::endl;
+   mf->Assemble();
+   mf->Finalize();
+   std::cout<<"M ready"<<std::endl;
+   cf->Assemble();
+   cf->Finalize();
+   std::cout<<"C ready"<<std::endl;
+   wf->Assemble();
+   wf->Finalize();
+   std::cout<<"W ready"<<std::endl;
+
+   std::cout<<"Cons"<<std::endl;
+
+   //(-M*\omega^2+i*\omega*C+K)u=f
+
+   unique_ptr<HypreParMatrix> kmat; kmat.reset(kf->ParallelAssemble());
+   unique_ptr<HypreParMatrix> cmat; cmat.reset(cf->ParallelAssemble());
+   unique_ptr<HypreParMatrix> mmat; mmat.reset(mf->ParallelAssemble());
+   unique_ptr<HypreParMatrix> wmat; wmat.reset(wf->ParallelAssemble());
+
+   kmat->EliminateBC(ess_dofs,Operator::DiagonalPolicy::DIAG_ONE);
+   mmat->EliminateBC(ess_dofs,Operator::DiagonalPolicy::DIAG_ZERO);
+   cmat->EliminateBC(ess_dofs,Operator::DiagonalPolicy::DIAG_ONE);
+   wmat->EliminateBC(ess_dofs,Operator::DiagonalPolicy::DIAG_ONE);
+
+   //kmat->EliminateRowsCols(ess_dofs); kmat->EliminateZeroRows();
+   //mmat->EliminateRowsCols(ess_dofs);
+   //cmat->EliminateRowsCols(ess_dofs);
+   //wmat->EliminateRowsCols(ess_dofs); wmat->EliminateZeroRows();
+
+
+   const Operator *P = kf->GetProlongation();
+   ConstrainedOperator ckf(new RAPOperator(*P,*kf,*P),ess_dofs,true, Operator::DiagonalPolicy::DIAG_ONE);
+   ConstrainedOperator cmf(new RAPOperator(*P,*mf,*P),ess_dofs,true, Operator::DiagonalPolicy::DIAG_ZERO);
+   ConstrainedOperator ccf(new RAPOperator(*P,*cf,*P),ess_dofs,true, Operator::DiagonalPolicy::DIAG_ONE);
+   ConstrainedOperator cwf(new RAPOperator(*P,*wf,*P),ess_dofs,true, Operator::DiagonalPolicy::DIAG_ONE);
+
+   OperatorHandle hkf;
+   OperatorHandle hmf;
+   OperatorHandle hcf;
+   OperatorHandle hwf;
+   kf->FormSystemMatrix(ess_dofs, hkf);
+   mf->FormSystemMatrix(ess_dofs, hmf);
+   cf->FormSystemMatrix(ess_dofs, hcf);
+   wf->FormSystemMatrix(ess_dofs, hwf);
+
+
+   unique_ptr<HypreParMatrix> bm;
+   {
+       Array2D<const HypreParMatrix*> am(2,2);
+       am(0,0)=kmat.get();
+       am(0,1)=cmat.get();
+       am(1,0)=cmat.get();
+       am(1,1)=kmat.get();
+
+       Array2D<real_t> cm(2,2);
+       cm(0,0)=1.0; cm(0,1)=1;
+       cm(1,0)=1.0; cm(1,1)=-1;
+
+       bm.reset(HypreParMatrixFromBlocks(am,&cm));
+    }
+
+   std::cout<<"Forcing!"<<std::endl;
+
+   //force coefficient
+   Vector A(dim); A=0.0; A[1]=1.0; //firce amplitude
+   Vector cc(dim); cc[0]=3.9; cc[1]=0.5; if(dim==3){cc[2]=0.0;} //force possition
+   VectorForceCoeff fc(A,cc,0.1);
+
+   ParLinearForm lf(vfes);
+   lf.AddDomainIntegrator(new VectorDomainLFIntegrator(fc));
+   lf.Assemble();
+
+   std::cout<<"LF is assmebled!"<<std::endl;
+
+   Array<int> block_true_offsets;
+   block_true_offsets.SetSize(3);
+   block_true_offsets[0] = 0;
+   block_true_offsets[1] = vfes->GetTrueVSize();
+   block_true_offsets[2] = vfes->GetTrueVSize();
+   block_true_offsets.PartialSum();
+
+   // solution
+   BlockVector x; x.Update(block_true_offsets); x=0.0;
+   // RHS
+   BlockVector f; f.Update(block_true_offsets);
+   std::cout<<pmesh.GetMyRank()<<" f0.size="<<f.GetBlock(0).Size()<<std::endl;
+   std::cout<<pmesh.GetMyRank()<<" f1.size="<<f.GetBlock(1).Size()<<std::endl;
+   std::cout<<pmesh.GetMyRank()<<" cwf.size="<<cwf.Width()<<std::endl;
+   std::cout<<pmesh.GetMyRank()<<" vfes.size="<<vfes->GetTrueVSize()<<std::endl;
+   lf.ParallelAssemble(f.GetBlock(0));
+   lf.ParallelAssemble(f.GetBlock(1));
+
+   std::cout<<"RHS is ready!"<<std::endl;
+
+   //cwf.EliminateRHS(x.GetBlock(0), f.GetBlock(0));
+   //cwf.EliminateRHS(x.GetBlock(1), f.GetBlock(1));
+
+   BlockOperator bop(block_true_offsets);
+   bop.SetBlock(0,0,wmat.get(),1.0);
+   bop.SetBlock(0,1,cmat.get(),1.0);
+   bop.SetBlock(1,0,cmat.get(),1.0);
+   bop.SetBlock(1,1,wmat.get(),-1.0);
+
+
+   SumOperator W(mmat.get(),1.0,mmat.get(),0.0,false,false);
+
+   ParLORDiscretization lork(*kf,ess_dofs);
+   HypreParMatrix& lorkm=lork.GetAssembledMatrix();
+
+   ParLORDiscretization lorc(*cf,ess_dofs);
+   HypreParMatrix& lorcm=lorc.GetAssembledMatrix();
+
+
+   ParFiniteElementSpace& lorfes=lork.GetParFESpace();
+
+   std::cout<<"Allocate PRESB"<<std::endl;
+   PRESBPrec* prec=new PRESBPrec(pmesh.GetComm(),1);
+   //prec->SetOperators(wmat.get(),mmat.get(),1.0,freq*freq,1, kmat.get());
+   //prec->SetOperators(kmat.get(),cmat.get(),1.0,1.0);
+   prec->SetOperators(&lorkm,&lorcm,1.0,1.0);
+   prec->SetAbsTol(1e-12);
+   prec->SetRelTol(1e-1);
+   prec->SetMaxIter(1000);
+
+
+   prec->Mult(f,x);
+   {
+       real_t rr=mfem::InnerProduct(pmesh.GetComm(), x,x);
+       if(pmesh.GetMyRank()==0){
+           std::cout<<"Residual="<<sqrt(rr)<<std::endl;
+       }
+   }
+
+
+   Vector xstat(f.GetBlock(0));xstat=0.0;
+   /*
+   {
+       //do the static solution
+       CGSolver ls(pmesh.GetComm());
+       HypreBoomerAMG amg(lorkm);
+       amg.SetElasticityOptions(vfes);
+       ls.SetOperator(*kmat);
+       ls.SetPreconditioner(amg);
+       ls.iterative_mode=false;
+       ls.SetAbsTol(1e-12);
+       ls.SetRelTol(1e-12);
+       ls.SetMaxIter(1000);
+       ls.SetPrintLevel(1);
+
+       ls.Mult(f.GetBlock(0),xstat);
+   }*/
+
+
+
+
+
+   ///eigenvalues check
+   {
+       CGSolver ls(pmesh.GetComm());
+       HypreBoomerAMG amg;
+       //amg.SetOperator(*kmat);
+       amg.SetOperator(lorkm);
+       //amg.SetElasticityOptions(&lorfes);
+       ls.SetOperator(*kmat);
+       ls.SetPreconditioner(amg);
+       ls.iterative_mode=false;
+       ls.SetPrintLevel(-1);
+       ls.SetAbsTol(1e-12);
+       ls.SetRelTol(1e-12);
+       ls.SetMaxIter(1000);
+
+       MUMPSSolver mumps(pmesh.GetComm());
+       mumps.SetPrintLevel(1);
+       mumps.SetMatrixSymType(MUMPSSolver::MatType::SYMMETRIC_POSITIVE_DEFINITE);
+       mumps.SetOperator(*kmat);
+
+       //ProductOperator pOp(&ls,mmat.get(),false,false);
+       LocProductOperator pOp(&mumps,mmat.get());
+
+       RandomizedSubspaceIteration ss(pmesh.GetComm());
+       ss.SetConstrDOFs(ess_dofs);
+       ss.SetNumModes(10);
+       ss.SetNumIter(10);
+       ss.SetOperator(pOp);
+       ss.Solve();
+
+       {
+
+           const std::vector<Vector>& vecs=ss.GetModes();
+           Vector rr(vecs[0]);
+
+
+           if(myrank==0){ std::cout<<std::endl;
+                          std::cout<<"Num modes="<<ss.GetNumModes()<<std::endl;}
+           for(int i=0;i<10;i++){
+                   kmat->Mult(vecs[i],rr);
+               for(int j=0;j<10;j++){
+                   real_t gp=InnerProduct (pmesh.GetComm(), vecs[j], rr);
+                   if(myrank==0){std::cout<<gp<<" ";}
+               }
+               if(myrank==0){std::cout<<std::endl;}
+           }
+
+           if(myrank==0){ std::cout<<std::endl;}
+           for(int i=0;i<10;i++){
+                   mmat->Mult(vecs[i],rr);
+               for(int j=0;j<10;j++){
+                   real_t gp=InnerProduct (pmesh.GetComm(), vecs[j], rr);
+                   if(myrank==0){std::cout<<gp<<" ";}
+               }
+               if(myrank==0){std::cout<<std::endl;}
+           }
+
+           if(myrank==0){ std::cout<<std::endl;}
+           for(int i=0;i<10;i++){
+                   kmat->Mult(vecs[i],rr);
+                   real_t gp=InnerProduct (pmesh.GetComm(), vecs[i], rr);
+
+                   mmat->Mult(vecs[i],rr);
+                   real_t kp=InnerProduct (pmesh.GetComm(), vecs[i], rr);
+
+                   if(myrank==0){std::cout<<"i="<<i<<" "<<kp/gp<<std::endl;}
+           }
+       }
+
+       AdaptiveRandomizedGenEig ae(pmesh.GetComm());
+       ae.SetOperators(*mmat,*kmat,mumps);
+       ae.SetNumModes(50);
+       ae.SetNumIter(4);
+       ae.SolveNA();
+
+       {
+
+           const std::vector<Vector>& vecs=ae.GetModes();
+           Vector rr(vecs[0]);
+
+
+           if(myrank==0){ std::cout<<std::endl;
+                          std::cout<<"Num modes="<<ss.GetNumModes()<<std::endl;}
+           for(int i=0;i<10;i++){
+                   kmat->Mult(vecs[i],rr);
+               for(int j=0;j<10;j++){
+                   real_t gp=InnerProduct (pmesh.GetComm(), vecs[j], rr);
+                   if(myrank==0){std::cout<<gp<<" ";}
+               }
+               if(myrank==0){std::cout<<std::endl;}
+           }
+
+           if(myrank==0){ std::cout<<std::endl;}
+           for(int i=0;i<10;i++){
+                   mmat->Mult(vecs[i],rr);
+               for(int j=0;j<10;j++){
+                   real_t gp=InnerProduct (pmesh.GetComm(), vecs[j], rr);
+                   if(myrank==0){std::cout<<gp<<" ";}
+               }
+               if(myrank==0){std::cout<<std::endl;}
+           }
+
+           if(myrank==0){ std::cout<<std::endl;}
+           for(int i=0;i<10;i++){
+                   mmat->Mult(vecs[i],rr);
+                   real_t gp=InnerProduct (pmesh.GetComm(), vecs[i], rr);
+
+                   kmat->Mult(vecs[i],rr);
+                   real_t kp=InnerProduct (pmesh.GetComm(), vecs[i], rr);
+
+                   //if(myrank==0){std::cout<<"i="<<i<<" "<<kp/gp<<std::endl;}
+                   if(myrank==0){std::cout<<kp/gp<<std::endl;}
+           }
+       }
+
+       {
+
+           mfem::ParGridFunction disp; disp.SetSpace(vfes);
+           const std::vector<Vector>& vecs=ae.GetModes();
+           ParaViewDataCollection paraview_dc("eigt", &pmesh);
+           paraview_dc.SetPrefixPath("ParaView");
+           paraview_dc.SetLevelsOfDetail(order);
+           paraview_dc.SetDataFormat(VTKFormat::BINARY);
+           paraview_dc.SetHighOrderOutput(true);
+           paraview_dc.RegisterField("disp",&disp);
+
+           for(int i=0;i<ae.GetNumModes();i++){
+               disp.SetFromTrueDofs(vecs[i]);
+               paraview_dc.SetCycle(i);
+               paraview_dc.SetTime(double(i));
+               paraview_dc.Save();
+           }
+       }
+
+
+
+
+/*
+       std::random_device rd;
+       std::mt19937 generator(rd());
+       // Create a normal distribution object
+       std::normal_distribution<real_t> distribution(0.0, 1.0);
+
+       std::vector<Vector> pvecs; pvecs.resize(5);
+       std::vector<Vector> ovecs; ovecs.resize(5);
+       for(int i=0;i<5;i++){
+           pvecs[i].SetSize(vecs[0].Size());
+           for(int j=0;j<vecs[0].Size();j++){
+               (pvecs[i])[j]=distribution(generator);
+           }
+       }
+
+
+       if(myrank==0){ std::cout<<std::endl;
+                      std::cout<<"Num modes="<<ss.GetNumModes()<<std::endl;}
+       for(int i=0;i<ss.GetNumModes();i++){
+               kmat->Mult(vecs[i],rr);
+           for(int j=0;j<ss.GetNumModes();j++){
+               real_t gp=InnerProduct (pmesh.GetComm(), vecs[j], rr);
+               if(myrank==0){std::cout<<gp<<" ";}
+           }
+           if(myrank==0){std::cout<<std::endl;}
+       }
+
+
+
+       ae.Ortho(kmat.get(), pvecs, ovecs);
+
+       if(myrank==0){ std::cout<<std::endl;}
+       for(int i=0;i<5;i++){
+           kmat->Mult(ovecs[i],rr);
+           for(int j=0;j<5;j++){
+               real_t gp=InnerProduct (pmesh.GetComm(), ovecs[j], rr);
+               if(myrank==0){std::cout<<gp<<" ";}
+           }
+           if(myrank==0){std::cout<<std::endl;}
+       }
+
+       if(myrank==0){std::cout<<std::endl;}
+       for(int i=0;i<5;i++){
+           mmat->Mult(pvecs[i],rr);
+           for(int j=0;j<5;j++){
+               real_t gp=InnerProduct (pmesh.GetComm(), pvecs[j], rr);
+               if(myrank==0){std::cout<<gp<<" ";}
+           }
+           if(myrank==0){std::cout<<std::endl;}
+       }
+
+       if(myrank==0){std::cout<<std::endl;}
+       ae.Ortho(mmat.get(),pvecs);
+       for(int i=0;i<5;i++){
+           mmat->Mult(pvecs[i],rr);
+           for(int j=0;j<5;j++){
+               real_t gp=InnerProduct (pmesh.GetComm(), pvecs[j], rr);
+               if(myrank==0){std::cout<<gp<<" ";}
+           }
+           if(myrank==0){std::cout<<std::endl;}
+       }
+       */
+
+   }
+
+
+   delete prec;
+   delete vfes;
+   delete vfec;
+   Mpi::Finalize();
+   return 0;
+
+   //set the linear solver
+   FGMRESSolver* ls=new FGMRESSolver(pmesh.GetComm());
+   //ls->SetOperator(bop);
+   ls->SetOperator(*bm);
+   ls->SetAbsTol(1e-12);
+   ls->SetRelTol(1e-12);
+   ls->SetMaxIter(1000);
+   ls->SetPrintLevel(1);
+   ls->SetKDim(100);
+
+   ls->SetOperator(bop);
+   ls->SetPreconditioner(*prec);
+   ls->Mult(f,x);
+   delete ls;
+
+   delete prec;
+
+   //check the solution
+   {
+       Vector xm(x); xm=0.0;
+       MUMPSSolver mumps(bm->GetComm());
+       mumps.SetPrintLevel(2);
+       mumps.SetMatrixSymType(MUMPSSolver::MatType::UNSYMMETRIC);
+       mumps.SetOperator(*bm);
+       mumps.Mult(f, xm);
+
+       xm.Add(-1,x);
+       real_t rr=mfem::InnerProduct(pmesh.GetComm(), xm,xm);
+       if(pmesh.GetMyRank()==0){
+           std::cout<<"|xm-x|="<<sqrt(rr)<<std::endl;
+       }
+   }
+
+
+   //check the solutions
+   ParGridFunction rx(vfes); rx.SetFromTrueDofs(x.GetBlock(0));
+   ParGridFunction ix(vfes); ix.SetFromTrueDofs(x.GetBlock(1)); ix*=-1.0;
+   ParGridFunction xs(vfes); xs.SetFromTrueDofs(xstat);
+   {
+        ParaViewDataCollection paraview_dc("stokes", &pmesh);
+        paraview_dc.SetPrefixPath("ParaView");
+        paraview_dc.SetLevelsOfDetail(order);
+        paraview_dc.SetDataFormat(VTKFormat::BINARY);
+        paraview_dc.SetHighOrderOutput(true);
+        paraview_dc.SetCycle(0);
+        paraview_dc.SetTime(0.0);
+        paraview_dc.RegisterField("re",&rx);
+        paraview_dc.RegisterField("im",&ix);
+        paraview_dc.RegisterField("xs",&xs);
+        paraview_dc.Save();
+   }
+
+   BlockVector r(f);
+
+   bop.Mult(x,r);
+   r.Add(-1,f);
+   real_t rr=mfem::InnerProduct(pmesh.GetComm(), r,r);
+   if(pmesh.GetMyRank()==0){
+       std::cout<<"Residual="<<sqrt(rr)<<std::endl;
+   }
+
+   delete vfes;
+   delete vfec;
+   Mpi::Finalize();
+   return 0;
+}

miniapps/mtop/test_stokes.cpp

@@ -0,0 +1,176 @@
+#include "mfem.hpp"
+#include "stokes_solver.hpp"
+
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+
+class GyroCoeff:public Coefficient
+{
+public:
+    GyroCoeff():l(0.0),scale(1.0){}
+
+    GyroCoeff(real_t l_,real_t s=1.0, real_t h=0.0)
+    {
+        l=l_;
+        scale=s;
+        shift=h;
+    }
+
+    virtual real_t Eval(ElementTransformation &T,
+                    const IntegrationPoint &ip)
+    {
+        real_t x[3];
+        Vector transip(x, 3);
+        T.Transform(ip, transip);
+
+        //real_t r=sin(l*x[0])*cos(l*x[1])+sin(l*x[1])*cos(l*x[2])+sin(l*x[2])*cos(l*x[0])-shift;
+        real_t r=cos(l*x[0])*cos(l*x[1])*cos(l*x[2])-shift;
+
+        if(r>=real_t(0.0))
+        {
+            if(x[0]<0.05){return 0.0;}
+            if(x[0]>1.45){return 0.0;}
+            if(x[1]<0.05){return 0.0;}
+            if(x[1]>0.95){return 0.0;}
+            return r*scale;
+        }else{
+            return 0.0;
+        }
+    }
+private:
+    real_t l;
+    real_t scale;
+    real_t shift;
+};
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI and HYPRE.
+   Mpi::Init();
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../../data/star.mesh";
+   int order = 1;
+   bool static_cond = false;
+   bool pa = false;
+   bool fa = false;
+   const char *device_config = "cpu";
+   bool visualization = true;
+   bool algebraic_ceed = false;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&fa, "-fa", "--full-assembly", "-no-fa",
+                  "--no-full-assembly", "Enable Full Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+#ifdef MFEM_USE_CEED
+   args.AddOption(&algebraic_ceed, "-a", "--algebraic",
+                  "-no-a", "--no-algebraic",
+                  "Use algebraic Ceed solver");
+#endif
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 3. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   if (myid == 0) { device.Print(); }
+
+   // 4. Read the (serial) mesh from the given mesh file on all processors.  We
+   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   //    and volume meshes with the same code.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   // 5. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement. We choose
+   //    'ref_levels' to be the largest number that gives a final mesh with no
+   //    more than 10,000 elements.
+   {
+      int ref_levels =
+         (int)floor(log(1000./mesh.GetNE())/log(2.)/dim);
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh.UniformRefinement();
+      }
+   }
+
+   // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+   {
+      int par_ref_levels = 1;
+      for (int l = 0; l < par_ref_levels; l++)
+      {
+         pmesh.UniformRefinement();
+      }
+   }
+   if(Mpi::WorldRank()==0){
+       std::cout<<pmesh.GetNE()<<std::endl;
+   }
+
+   GyroCoeff gc(6.0*M_PI,1000.0,0.3);
+
+   StokesOperator* sp=new StokesOperator(&pmesh,2);
+   sp->AddVelocityBC(1,4,0.0);
+   sp->AddVelocityBC(2,0,1.0);
+
+   mfem::ConstantCoefficient mu(0.01);
+   sp->SetViscosoty(mu);
+   sp->SetBrinkman(gc);
+   sp->Assemble();
+
+   sp->FSolve();
+
+   {
+       ParGridFunction brink(sp->GetPressure()); brink.ProjectCoefficient(gc);
+       ParaViewDataCollection paraview_dc("stokes", &pmesh);
+       paraview_dc.SetPrefixPath("ParaView");
+       paraview_dc.SetLevelsOfDetail(order);
+       paraview_dc.SetDataFormat(VTKFormat::BINARY);
+       paraview_dc.SetHighOrderOutput(true);
+       paraview_dc.SetCycle(0);
+       paraview_dc.SetTime(0.0);
+       paraview_dc.RegisterField("velo",&(sp->GetVelocity()));
+       paraview_dc.RegisterField("pres",&(sp->GetPressure()));
+       paraview_dc.RegisterField("brink",&(brink));
+       paraview_dc.Save();
+   }
+
+   delete sp;
+   Mpi::Finalize();
+   return 0;
+}
+