mod

miniapps/mtop/CMakeLists.txt

@@ -9,46 +9,31 @@
 # terms of the BSD-3 license. We welcome feedback and contributions, see file
 # 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)
-
-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)
+        mtop_solvers.cpp
+        grain_reader.cpp)
 list(APPEND PARMTOP_COMMON_HEADERS
-	pparamnonlinearform.hpp)
+        mtop_solvers.hpp
+        grain_reader.hpp)
 
 convert_filenames_to_full_paths(PARMTOP_COMMON_SOURCES)
 convert_filenames_to_full_paths(PARMTOP_COMMON_HEADERS)
 
 set(PARMTOP_COMMON_FILES
-	EXTRA_SOURCES ${PARMTOP_COMMON_SOURCES} ${SEQMTOP_COMMON_SOURCES}
-	EXTRA_HEADERS ${PARMTOP_COMMON_HEADERS} ${SEQMTOP_COMMON_HEADERS})
+        EXTRA_SOURCES ${PARMTOP_COMMON_SOURCES} ${SEQMTOP_COMMON_SOURCES}
+        EXTRA_HEADERS ${PARMTOP_COMMON_HEADERS} ${SEQMTOP_COMMON_HEADERS})
 
-# message(STATUS "PARMTOP_COMMON_FILES: ${PARMTOP_COMMON_FILES}")
-# message(STATUS "SEQMTOP_COMMON_FILES: ${SEQMTOP_COMMON_FILES}")
 
-add_mfem_miniapp(parheat
-     MAIN parheat.cpp
-     ${PARMTOP_COMMON_FILES}
-     LIBRARIES mfem)
+add_mfem_miniapp(test_stokes
+       MAIN test_stokes.cpp
+       ${PARMTOP_COMMON_FILES}
+       LIBRARIES mfem)
+
+add_mfem_miniapp(test_grain
+       MAIN test_grain.cpp
+       ${PARMTOP_COMMON_FILES}
+       LIBRARIES mfem)
 
 endif ()

miniapps/mtop/grain_reader.cpp

@@ -0,0 +1,419 @@
+// Copyright (c) 2010-2021, 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 <fstream>
+#include <iostream>
+
+#include "grain_reader.hpp"
+#include "mfem.hpp"
+
+
+namespace mfem
+{
+
+    //------------------------------------------------------------------------------
+
+    bool
+    file_exists( const std::string & aPath )
+    {
+        // test if file exists
+        std::ifstream tFile( aPath );
+
+        // save result into output variable
+        bool aFileExists;
+
+        if( tFile )
+        {
+            // close file
+            tFile.close();
+            aFileExists = true;
+        }
+        else
+        {
+            aFileExists = false;
+        }
+
+        return aFileExists;
+    }
+
+    //------------------------------------------------------------------------------
+
+    std::vector<std::string> split_string(
+            const std::string  & aString,
+            const std::string  & aDelim)
+    {
+        // create empty cell of strings
+        std::vector<std::string> VectorOfStrings;
+
+        size_t start;
+        size_t end = 0;
+
+        while ((start = aString.find_first_not_of(aDelim, end)) != std::string::npos)
+        {
+            end = aString.find(aDelim, start);
+            VectorOfStrings.push_back(aString.substr(start, end - start));
+        }
+
+        return VectorOfStrings;
+    }
+
+    //------------------------------------------------------------------------------
+
+    Ascii::Ascii( const std::string & aPath, const FileMode & aMode ) :
+                mMode( aMode )
+    {
+        // test if path is absolute
+        if( aPath.substr( 0,1 ) == "/" )
+        {
+            mPath = aPath;
+        }
+        // test if path is relative
+        else if( aPath.substr( 0,2 ) == "./" )
+        {
+            mPath = aPath;
+        }
+        else
+        {
+            MFEM_ASSERT( false, "");
+            //mPath = std::sprint( "%s/%s", std::getenv( "PWD" ), aPath.c_str() );
+        }
+
+        switch ( aMode )
+        {
+            case( FileMode::OPEN_RDONLY ) :
+            {
+                this->load_buffer();
+                break;
+            }
+            case( FileMode::NEW ) :
+            {
+                mBuffer.clear();
+                break;
+            }
+            default:
+            {
+                MFEM_ASSERT( false, "Unknown file mode for ASCII file" );
+            }
+        }
+    }
+
+    //------------------------------------------------------------------------------
+
+    Ascii::~Ascii()
+    {
+        MFEM_ASSERT( ! mChangedSinceLastSave, "File  was changed but never saved." );
+
+        mBuffer.clear();
+    }
+    //------------------------------------------------------------------------------
+
+    bool Ascii::save()
+    {
+        MFEM_ASSERT( mMode != FileMode::OPEN_RDONLY,
+                "File can't be saved since it is opened in write protected mode." );
+
+        // open file
+        std::ofstream tFile( mPath.c_str(),  std::ofstream::trunc );
+
+        if( tFile )
+        {
+            // save buffer to file
+            for( std::string & tLine : mBuffer )
+            {
+                tFile << tLine << std::endl;
+            }
+
+            tFile.close();
+        }
+        else
+        {
+            MFEM_ASSERT( false,
+                    "Something went wrong while trying to save." );
+        }
+
+        mChangedSinceLastSave = false;
+
+        return mChangedSinceLastSave;
+    }
+
+    //------------------------------------------------------------------------------
+
+    int Ascii::length() const
+    {
+        return mBuffer.size();
+    }
+
+    //------------------------------------------------------------------------------
+
+    std::string & Ascii::line( const int aLineNumber )
+    {
+        return mBuffer[aLineNumber];
+    }
+
+    //------------------------------------------------------------------------------
+
+    const std::string & Ascii::line( const int aLineNumber ) const
+    {
+        return mBuffer[aLineNumber];
+    }
+
+    //------------------------------------------------------------------------------
+
+    void Ascii::print( const std::string & aLine )
+    {
+        mBuffer.push_back( aLine );
+
+        mChangedSinceLastSave = true;
+    }
+
+    //------------------------------------------------------------------------------
+
+    void Ascii::load_buffer()
+    {
+        // tidy up buffer
+        mBuffer.clear();
+
+        // make sure that file exists
+        MFEM_ASSERT( file_exists( mPath ),
+                "File does not exist." );
+
+        // open file
+        std::ifstream tFile( mPath );
+
+        // test if file can be opened
+        if( tFile )
+        {
+            // temporary container for string
+            std::string tLine;
+
+            while ( std::getline( tFile, tLine ) )
+            {
+                mBuffer.push_back( tLine );
+            }
+
+            // close file
+            tFile.close();
+        }
+        else
+        {
+            MFEM_ASSERT( false, "Someting went wrong while opening file\n " );
+        }
+    }
+
+    //------------------------------------------------------------------------------
+
+    GrainReader::GrainReader( mfem::ParMesh * mesh, std::string & name ) 
+      : mesh_(mesh), name_(name)
+   {
+      dim = mesh->Dimension();
+
+      int tNumVertices  = mesh->GetNV();
+      for (int i = 0; i < tNumVertices; ++i)
+      {
+         double * Coords = mesh->GetVertex(i);
+
+         xMax = std::max(xMax, Coords[ 0 ]);
+         yMax = std::max(yMax, Coords[ 1 ]);
+         zMax = std::max(zMax, Coords[ 2 ]);
+
+         xMin = std::min(xMin, Coords[ 0 ]);
+         yMin = std::min(yMin, Coords[ 1 ]);
+         zMin = std::min(zMin, Coords[ 2 ]);
+      }
+
+      MPI_Allreduce(MPI_IN_PLACE, &xMax, 1, MPI_DOUBLE, MPI_MAX, mesh_->GetComm());
+      MPI_Allreduce(MPI_IN_PLACE, &yMax, 1, MPI_DOUBLE, MPI_MAX, mesh_->GetComm());
+      MPI_Allreduce(MPI_IN_PLACE, &zMax, 1, MPI_DOUBLE, MPI_MAX, mesh_->GetComm());
+
+      MPI_Allreduce(MPI_IN_PLACE, &xMin, 1, MPI_DOUBLE, MPI_MIN, mesh_->GetComm());
+      MPI_Allreduce(MPI_IN_PLACE, &yMin, 1, MPI_DOUBLE, MPI_MIN, mesh_->GetComm());
+      MPI_Allreduce(MPI_IN_PLACE, &zMin, 1, MPI_DOUBLE, MPI_MIN, mesh_->GetComm());
+
+      Lx = xMax - xMin;
+      Ly = yMax - yMin;
+      Lz = zMax - zMin;
+
+      this->readGrainFile();
+   }
+
+   void GrainReader::readGrainFile()
+   {
+   mfem::Ascii tAsciiReader( name_, FileMode::OPEN_RDONLY );
+
+   int tNumLines = tAsciiReader.length();
+
+   particalType.reserve(9*tNumLines);
+   xPos.reserve(9*tNumLines);
+   yPos.reserve(9*tNumLines);
+   zPos.reserve(9*tNumLines);
+   rad.reserve(9*tNumLines);
+
+   real_t maxRad = 0.0;
+
+   for( int Ik = 0; Ik < tNumLines; Ik++ )
+   {
+      const std::string & tFileLine = tAsciiReader.line( Ik );
+
+      std::vector<std::string> ListOfStrings = split_string( tFileLine, " " );
+
+      particalType.push_back(std::stod( ListOfStrings[1] ));
+      xPos        .push_back(std::stod( ListOfStrings[2] ));
+      yPos        .push_back(std::stod( ListOfStrings[3] ));
+      zPos        .push_back(std::stod( ListOfStrings[4] ));
+      rad         .push_back(std::stod( ListOfStrings[5] ) / 2.0);
+
+      maxRad = std::max( rad[Ik], maxRad );
+   }
+
+   maxRad += 1e-6;
+
+   for( int Ik = 0; Ik < tNumLines; Ik++ )
+   {
+      int pType = particalType[Ik];
+      real_t xCopy = xPos[Ik];
+      real_t yCopy = yPos[Ik];
+      real_t zCopy = zPos[Ik];
+      real_t radCopy = rad[Ik];
+
+      real_t xC;
+      real_t yC;
+
+      bool isCopyX = false;
+      bool isCopyY = false;
+      bool isCopyCorner = false;
+
+      if(xPos[Ik]  < (xMin + maxRad)){ 
+         xC = xCopy + Lx; 
+         isCopyX = true; }
+      else if(xPos[Ik]  > (xMin + Lx - maxRad)) { 
+         xC = xCopy - Lx; 
+         isCopyX = true; }
+
+      if(isCopyX)
+      {
+         particalType.push_back(pType);
+         xPos        .push_back(xC);
+         yPos        .push_back(yCopy);
+         zPos        .push_back(zCopy);
+         rad         .push_back(radCopy);
+      }
+
+      if(yPos[Ik] < (yMin + maxRad))      { 
+         yC = yCopy + Ly;
+         isCopyY = true; }
+      else if(yPos[Ik]  > (yMin + Ly - maxRad)) {
+         yC = yCopy - Ly;
+         isCopyY = true; }
+
+      if(isCopyY)
+      {
+         particalType.push_back(pType);
+         xPos        .push_back(xCopy);
+         yPos        .push_back(yC);
+         zPos        .push_back(zCopy);
+         rad         .push_back(radCopy);
+      }
+
+      if(xPos[Ik]  < (xMin + maxRad) && yPos[Ik]  < (yMin + maxRad))     { 
+         xC = xCopy + Lx;  
+         yC = yCopy + Ly;
+         isCopyCorner = true; }
+      else if(xPos[Ik]  < (xMin + maxRad) && yPos[Ik]  > (yMin + Ly - maxRad)) { 
+         xC = xCopy + Lx;    
+         yC = yCopy - Ly;
+         isCopyCorner = true; }
+
+      else if(xPos[Ik]  > (xMin + Lx - maxRad) && yPos[Ik]  < (yMin + maxRad))     { 
+         xC = xCopy - Lx;  
+         yC = yCopy + Ly;
+         isCopyCorner = true; }
+      else if(xPos[Ik]  > (xMin + Lx - maxRad) && yPos[Ik]  > (yMin + Ly - maxRad)) { 
+         xC = xCopy - Lx;  
+         yC = yCopy - Ly;
+         isCopyCorner = true; }
+
+      if(isCopyCorner)
+      {
+         particalType.push_back(pType);
+         xPos        .push_back(xC);
+         yPos        .push_back(yC);
+         zPos        .push_back(zCopy);
+         rad         .push_back(radCopy);
+      }
+   }
+
+   particalType.shrink_to_fit();
+   xPos        .shrink_to_fit();
+   yPos        .shrink_to_fit();
+   zPos        .shrink_to_fit();
+   rad         .shrink_to_fit();
+
+   numParticles = particalType.size();
+   }
+
+   void GrainReader::computeGridFunction( ::mfem::ParFiniteElementSpace& feSpace)
+   {
+      grainLSField.SetSpace(&feSpace);
+
+      int numNodes   = grainLSField.Size();
+      mfem::Vector locationVector(dim);
+
+      int numEle = mesh_->GetNE();
+
+      for ( int e = 0; e<numEle; e++)
+      {
+         const IntegrationRule &ir = feSpace.GetFE(e)->GetNodes();
+
+         // Transformation of the element with the pos_mesh coordinates.
+         mfem::IsoparametricTransformation Tr;
+         feSpace.GetElementTransformation(e, &Tr);
+
+         mfem::DenseMatrix pos_nodes;
+         Tr.Transform(ir, pos_nodes);
+         mfem::Vector valVec(pos_nodes.NumCols());
+
+         for ( int Ik = 0; Ik< pos_nodes.NumCols(); Ik++)
+         {
+            double LSVal = -1000.0;
+            for (int ii = 0; ii < numParticles; ii++)
+            {
+               double val = rad[ii] - pow(pow(std::abs(pos_nodes(0,Ik) - xPos[ii]), 2) 
+                  + pow(std::abs(pos_nodes(1,Ik) - yPos[ii]), 2) + pow(std::abs(pos_nodes(2,Ik) - zPos[ii]), 2), 0.5);
+
+               LSVal = std::max(val, LSVal);
+            }
+            valVec[Ik]= LSVal;
+         }
+
+         mfem::Array< int > dofs;
+         feSpace.GetElementDofs( e, dofs );
+
+         grainLSField.SetSubVector(dofs, valVec);
+      }
+
+      // for ( int Ik = 0; Ik<numNodes; Ik++)
+      // {
+      //    mesh_->GetNode(Ik, &locationVector[0]);
+      //    const double * pCoords(locationVector.GetData());
+
+      //    double LSVal = -1000.0;
+      //    for (int ii = 0; ii < numParticles; ii++)
+      //    {
+      //       double val = rad[ii] - pow(pow(std::abs(pCoords[0] - xPos[ii]), 2) 
+      //            + pow(std::abs(pCoords[1] - yPos[ii]), 2) + pow(std::abs(pCoords[2] - zPos[ii]), 2), 0.5);
+
+      //       LSVal = std::max(val, LSVal);
+      //    }
+      //    grainLSField[Ik]= LSVal;
+      // }
+   }
+
+}

miniapps/mtop/grain_reader.hpp

@@ -0,0 +1,199 @@
+// Copyright (c) 2010-2021, 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 PROJECTS_ASCII_HPP_
+#define PROJECTS_ASCII_HPP_
+
+#include <vector>
+#include <iomanip>
+#include <iostream>
+#include <sstream>
+
+#include <fstream>
+#include "mfem.hpp"
+
+namespace mfem
+{
+    //------------------------------------------------------------------------------
+
+    enum class FileMode
+    {
+            NEW,
+            OPEN_RDONLY,
+            OPEN_RDWR
+    };
+
+    //------------------------------------------------------------------------------
+
+    /**
+     * this function tests if a file exists
+     * @param aPath
+     * @return
+     */
+    bool file_exists( const std::string & aPath );
+
+
+    std::vector<std::string> split_string(
+            const std::string  & aString,
+            const std::string  & aDelim);
+
+    class Ascii
+    {
+            //------------------------------------------------------------------------------
+        protected:
+            //------------------------------------------------------------------------------
+
+            std::string          mPath;
+            const FileMode       mMode;
+            std::vector< std::string >       mBuffer;
+
+            bool                 mChangedSinceLastSave = false;
+
+            //------------------------------------------------------------------------------
+        public:
+            //------------------------------------------------------------------------------
+
+            Ascii(
+                    const std::string   & aPath,
+                    const enum FileMode & aMode );
+
+            //------------------------------------------------------------------------------
+
+            ~Ascii();
+
+            //------------------------------------------------------------------------------
+
+            // save the buffer to the file
+            bool save();
+
+            //------------------------------------------------------------------------------
+
+            void print( const std::string & aLine );
+
+            //------------------------------------------------------------------------------
+
+            /**
+             * return the number of lines
+             */
+            int length() const;
+
+            //------------------------------------------------------------------------------
+
+            std::string & line( const int aLineNumber );
+
+            //------------------------------------------------------------------------------
+
+            const std::string & line( const int aLineNumber ) const;
+
+            //------------------------------------------------------------------------------
+        private:
+            //------------------------------------------------------------------------------
+
+            void load_buffer();
+
+            //------------------------------------------------------------------------------
+    };
+
+    //------------------------------------------------------------------------------
+
+        template<typename T>
+        inline std::string stringify(T aValue)
+        {
+            std::ostringstream out;
+            out << aValue;
+            return out.str();
+        }
+
+        template<>
+        inline std::string stringify<bool>(bool aValue)
+        {
+            std::ostringstream out;
+            out << std::boolalpha << aValue;
+            return out.str();
+        }
+
+        template<>
+        inline std::string stringify<double>(double aValue)
+        {
+            std::ostringstream out;
+            out << std::setprecision(14) << std::scientific << aValue;
+            return out.str();
+        }
+
+        template<>
+        inline std::string stringify<long double>(long double aValue)
+        {
+            std::ostringstream out;
+            out << std::setprecision(14) << std::scientific << aValue;
+            return out.str();
+        }
+
+        template<>
+        inline std::string stringify<float>(float aValue)
+        {
+            std::ostringstream out;
+            out << std::setprecision(14) << std::scientific << aValue;
+            return out.str();
+        }
+
+
+class GrainReader
+{
+private:
+   mfem::ParMesh * mesh_;
+   std::string name_;
+   int dim;
+
+   double xMax = DBL_MIN;
+   double yMax = DBL_MIN;
+   double zMax = DBL_MIN;
+   double xMin = DBL_MAX;
+   double yMin = DBL_MAX;
+   double zMin = DBL_MAX;
+
+   double Lx, Ly, Lz;
+
+   std::vector<int> particalType; 
+   std::vector<real_t> xPos;
+   std::vector<real_t> yPos;
+   std::vector<real_t> zPos;
+   std::vector<real_t> rad;
+
+   ::mfem::ParGridFunction grainLSField;
+   ::mfem::QuadratureFunction grainQFField;
+
+   int numParticles;
+
+public:
+   GrainReader( mfem::ParMesh * mesh, std::string & name );
+
+   ~GrainReader(){}
+
+   void readGrainFile();
+
+   void computeGridFunction( ::mfem::ParFiniteElementSpace& feSpace);
+
+   void computeQuadratureFunction( QuadratureSpaceBase &qspace)
+   {
+      ::mfem::mfem_error("not implemented yet");
+      grainQFField.SetSpace(&qspace);
+   }
+
+   const ::mfem::ParGridFunction & getGrainGridFunction() const
+   { return grainLSField; }
+
+   const ::mfem::QuadratureFunction & getGrainQuadratureFunction() const
+   { return grainQFField; }
+};
+
+}
+
+#endif /* PROJECTS_ASCII_HPP_ */

miniapps/mtop/makefile

@@ -24,12 +24,12 @@ include $(DEFAULTS_MK)
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
-MTOP_COMMON_SRC = mtop_integrators.cpp paramnonlinearform.cpp pparamnonlinearform.cpp
+MTOP_COMMON_SRC = mtop_solvers.cpp
 
 MTOP_COMMON_OBJ = $(MTOP_COMMON_SRC:.cpp=.o)
 
-SEQ_MINIAPPS = seqheat
-PAR_MINIAPPS = parheat
+SEQ_MINIAPPS =
+PAR_MINIAPPS = test_stokes
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
 else

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,339 @@
+// 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 <memory>
+
+#include "mtop_solvers.hpp"
+
+using namespace mfem;
+
+StokesSolver::StokesSolver(ParMesh* mesh, int order_, bool zero_mean_press_):
+   pmesh(mesh),
+   order(order_),
+   dim(mesh->SpaceDimension()),
+   zero_mean_press(zero_mean_press_)
+{
+   if (order_<2) { order=2;}
+
+   vfec=new H1_FECollection(order, pmesh->Dimension());
+   pfec=new H1_FECollection(order-1, pmesh->Dimension());
+   vfes=new ParFiniteElementSpace(pmesh, vfec, pmesh->Dimension());
+   pfes=new ParFiniteElementSpace(pmesh, pfec);
+
+   vel.SetSpace(vfes); vel=0.0;
+   pre.SetSpace(pfes); pre=0.0;
+
+   avel.SetSpace(vfes); avel=0.0;
+   apre.SetSpace(pfes); apre=0.0;
+
+   brink.reset();
+   visc.reset(new ConstantCoefficient(0.001));
+
+   onecoeff.constant = 1.0;
+   zerocoef.constant = 0.0;
+
+   siz_u=vfes->TrueVSize();
+   siz_p=pfes->TrueVSize();
+
+   block_true_offsets.SetSize(3);
+   block_true_offsets[0] = 0;
+   block_true_offsets[1] = siz_u;
+   block_true_offsets[2] = siz_p;
+   block_true_offsets.PartialSum();
+   //set the width and the height of the operator
+   this->width=  block_true_offsets[2];
+   this->height= block_true_offsets[2];
+
+
+   sol.Update(block_true_offsets); sol=0.0;
+   rhs.Update(block_true_offsets); rhs=0.0;
+   adj.Update(block_true_offsets); adj=0.0;
+
+   ess_tdofv.SetSize(0);
+
+   bf11.reset();
+   bf12.reset();
+   bf21.reset();
+
+   SetLinearSolver();
+}
+
+StokesSolver::~StokesSolver()
+{
+   delete pfes;
+   delete vfes;
+   delete pfec;
+   delete vfec;
+}
+
+void StokesSolver::SetEssTDofsV(mfem::Array<int>& ess_dofs)
+{
+   // Set the essential boundary conditions
+   ess_dofs.DeleteAll();
+
+   Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+   ess_bdr=0;
+   for (auto it=vel_bcs.begin(); it!=vel_bcs.end(); ++it)
+   {
+      int attr = it->first;
+      ess_bdr[attr-1] = 1;
+   }
+   vfes->GetEssentialTrueDofs(ess_bdr,ess_dofs);
+}
+
+void StokesSolver::SetEssTDofsV(Vector& v) const
+{
+   for (auto it=vel_bcs.begin(); it!=vel_bcs.end(); ++it)
+   {
+      int attr = it->first;
+      std::shared_ptr<VectorCoefficient> coeff = it->second;
+      coeff->SetTime(real_t(0.0));
+
+      Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+      ess_bdr=0;
+      ess_bdr[attr-1] = 1;
+
+      mfem::Array<int> loc_tdofs;
+      vfes->GetEssentialTrueDofs(ess_bdr,loc_tdofs);
+      vel.ProjectBdrCoefficient(*coeff,ess_bdr);
+      vel.SetTrueVector();
+
+      // copy values to v
+      Vector &tvel=vel.GetTrueVector();
+      //vel.GetTrueDofs(tvel);
+      for (int j=0; j<loc_tdofs.Size(); j++)
+      {
+         v[loc_tdofs[j]]=tvel[loc_tdofs[j]];
+      }
+   }
+}
+
+void StokesSolver::SetEssVBC(ParGridFunction& pgf)
+{
+
+   for (auto it=vel_bcs.begin(); it!=vel_bcs.end(); ++it)
+   {
+      int attr = it->first;
+      std::shared_ptr<VectorCoefficient> coeff = it->second;
+      coeff->SetTime(real_t(0.0));
+      Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+      ess_bdr=0;
+      ess_bdr[attr-1] = 1;
+
+      pgf.ProjectBdrCoefficient(*coeff,ess_bdr);
+   }
+}
+
+void StokesSolver::DeleteBC()
+{
+   vel_bcs.clear();
+   ess_tdofv.DeleteAll();
+   //delete allocated matrices and forms
+
+   bf11.reset();
+   bf12.reset();
+   bf21.reset();
+
+   A11.reset();
+   A12.reset();
+   A21.reset();
+
+   A11e.reset();
+   A12e.reset();
+   A21e.reset();
+
+   bop.reset();
+   ls.reset();
+   prec.reset();
+}
+
+void StokesSolver::Assemble()
+{
+   //set BC
+   vel=real_t(0.0);
+   pre=real_t(0.0);
+   SetEssVBC(vel);
+   SetEssTDofsV(ess_tdofv);
+
+   //assemble block 11
+   bf11.reset(new ParBilinearForm(vfes));
+   bf11->AddDomainIntegrator(new ElasticityIntegrator(zerocoef,*visc));
+   if (nullptr!=brink.get())
+   {
+      bf11->AddDomainIntegrator(new VectorMassIntegrator(*brink));
+   }
+   bf11->Assemble(0);
+   bf11->Finalize(0);
+   A11.reset(bf11->ParallelAssemble());
+
+   //assemble block 12
+   bf12.reset(new ParMixedBilinearForm(pfes, vfes));
+   //bf12->AddDomainIntegrator(new GradientIntegrator());
+   bf12->AddDomainIntegrator(
+               new TransposeIntegrator(
+                   new VectorDivergenceIntegrator()));
+   bf12->Assemble(0);
+   bf12->Finalize(0);
+   A12.reset(bf12->ParallelAssemble());
+
+   //assemble block 21
+   bf21.reset(new ParMixedBilinearForm(vfes, pfes));
+   bf21->AddDomainIntegrator(
+               new VectorDivergenceIntegrator());
+   bf21->Assemble(0);
+   bf21->Finalize(0);
+   A21.reset(bf21->ParallelAssemble());
+
+   //set BC to the operators
+   A11e.reset(A11->EliminateRowsCols(ess_tdofv));
+   A12->EliminateRows(ess_tdofv);
+   A21e.reset(A21->EliminateCols(ess_tdofv));
+
+   //set the block operator
+   bop.reset(new BlockOperator(block_true_offsets));
+   bop->SetBlock(0,0,A11.get());
+   bop->SetBlock(0,1,A12.get());
+   bop->SetBlock(1,0,A21.get());
+
+   if (zero_mean_press)
+   {
+      V.SetSize(pfes->GetTrueVSize()); V=0.0;
+      ParLinearForm lf(pfes);
+      lf.AddDomainIntegrator(new DomainLFIntegrator(onecoeff));
+      lf.Assemble();
+      lf.ParallelAssemble(V);
+   }
+
+   //set the solver to GMRES
+   {
+      //GMRESSolver* gmres=new GMRESSolver(pmesh->GetComm());
+
+      //MINRESSolver* gmres=new MINRESSolver(pmesh->GetComm());
+      FGMRESSolver* gmres=new FGMRESSolver(pmesh->GetComm());
+      gmres->SetKDim(100);
+      gmres->SetRelTol(linear_rtol);
+      gmres->SetAbsTol(linear_atol);
+      gmres->SetMaxIter(linear_iter);
+      gmres->SetOperator(*bop);
+      gmres->SetPrintLevel(1);
+
+      //prec.reset(new DLSCPrec(A11.get(),A21.get(),A12.get(), zero_mean_press));
+
+      LSCStokesPrec* lsc=new LSCStokesPrec(vfes,pfes,visc,brink,ess_tdofv,
+                                           A11.get(),A12.get(),A21.get(),zero_mean_press);
+
+      prec.reset(lsc);
+      prec->SetMaxIter(100);
+      prec->SetAbsTol(1e-12);
+      prec->SetRelTol(1e-5);
+      gmres->SetPreconditioner(*prec);
+
+      ls.reset(gmres);
+   }
+
+}
+
+void StokesSolver::FSolve()
+{
+   Vector& vsol=sol.GetBlock(0);
+   Vector& psol=sol.GetBlock(1);
+
+   Vector& vrhs=rhs.GetBlock(0);
+   Vector& prhs=rhs.GetBlock(1);
+
+   //assemble the RHS
+   rhs=0.0;
+   if (nullptr!=vol_force.get())
+   {
+      ParLinearForm lf(vfes);
+      lf.AddDomainIntegrator(new VectorDomainLFIntegrator(*vol_force));
+      lf.Assemble();
+      lf.ParallelAssemble(vrhs);
+   }
+
+   //set the velocity BCs
+   SetEssTDofsV(vsol);
+
+
+   //modify the rhs
+   A21e->Mult(-1.0,vsol,1.0,prhs);
+   A11->EliminateBC(*A11e,ess_tdofv,vsol,vrhs);
+
+   //solve the linear system
+   ls->Mult(rhs,sol);
+
+}
+
+void StokesSolver::ASolve(mfem::Vector &rhs)
+{
+   MultTranspose(rhs,adj);
+}
+
+void StokesSolver::Mult(const mfem::Vector &x, mfem::Vector &y) const
+{
+   //copy x to rhs
+   {
+      int N = x.Size();
+      const real_t *xp = x.Read();
+      real_t *rp = rhs.ReadWrite();
+      forall(N, [=] MFEM_HOST_DEVICE(int i) { rp[i] = xp[i]; });
+   }
+
+   BlockVector yb(y, block_true_offsets);
+   Vector& vsol=yb.GetBlock(0);
+   Vector& psol=yb.GetBlock(1);
+
+   Vector& vrhs=rhs.GetBlock(0);
+   Vector& prhs=rhs.GetBlock(1);
+
+   //set the velocity BCs
+   SetEssTDofsV(vsol);
+
+   //modify the rhs
+   A21e->Mult(-1.0,vsol,1.0,prhs);
+   A11->EliminateBC(*A11e,ess_tdofv,vsol,vrhs);
+
+   //solve the linear system
+   ls->Mult(rhs,yb);
+}
+
+void StokesSolver::MultTranspose(const mfem::Vector &x, mfem::Vector &y) const
+{
+   //copy x to rhs
+   {
+      int N = x.Size();
+      const real_t *xp = x.Read();
+      real_t *rp = rhs.ReadWrite();
+      forall(N, [=] MFEM_HOST_DEVICE(int i) { rp[i] = xp[i]; });
+   }
+
+   BlockVector yb(y, block_true_offsets);
+   Vector& vsol=yb.GetBlock(0);
+   Vector& psol=yb.GetBlock(1);
+
+   Vector& vrhs=rhs.GetBlock(0);
+   Vector& prhs=rhs.GetBlock(1);
+
+   //set zero velocity bc
+   {
+      int N = ess_tdofv.Size();
+      real_t *yp = vsol.ReadWrite();
+      const int *ep = ess_tdofv.Read();
+      forall(N, [=] MFEM_HOST_DEVICE(int i) { yp[ep[i]] = 0.0; });
+   }
+
+   //modify the rhs
+   //A21e->Mult(-1.0,vsol,1.0,prhs); // vsol is zero at the BC
+   A11->EliminateBC(*A11e,ess_tdofv,vsol,vrhs);
+
+   //solve the linear system
+   ls->Mult(rhs,yb);
+}

miniapps/mtop/mtop_solvers.hpp

@@ -0,0 +1,432 @@
+// 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.
+#pragma once
+
+#include <memory>
+
+#include "mfem.hpp"
+
+using real_t = mfem::real_t;
+
+class StokesSolver:public mfem::Operator
+{
+public:
+   StokesSolver(mfem::ParMesh* mesh, int order_, bool zero_mean_press_=false);
+
+   virtual
+   ~StokesSolver();
+
+   /// Set the Linear Solver
+   void SetLinearSolver(const real_t rtol = 1e-8,
+                        const real_t atol = 1e-12,
+                        const int miter = 200)
+   {
+      linear_atol=atol;
+      linear_rtol=rtol;
+      linear_iter=miter;
+   }
+
+   /// Sets BC dofs, bilinear form, preconditioner and solver.
+   /// Should be called before calling Mult of MultTranspose
+   virtual void Assemble();
+
+   /// Sets Brinkman coefficient
+   void SetBrink(std::shared_ptr<mfem::Coefficient> br_){
+       brink=br_;
+   }
+
+   /// Sets viscosity
+   void SetVisc(std::shared_ptr<mfem::Coefficient> vs_){
+       visc=vs_;
+   }
+
+   /// Solves the forward problem.
+   void FSolve();
+
+   /// Solves the adjoint with the provided rhs.
+   void ASolve(mfem::Vector &rhs);
+
+   /// Clear all  BC
+   void DeleteBC();
+
+   /// Set the values of the volumetric force.
+   void SetVolForce(real_t fx, real_t fy, real_t fz = 0.0);
+
+   //Set zero mean pressure BC
+   void SetZeroMeanPressure(bool val_=true)
+   {
+      zero_mean_press=val_;
+   }
+
+   //velocity boundary conditions
+   void AddVelocityBC(int id, std::shared_ptr<mfem::VectorCoefficient> val)
+   {
+      vel_bcs[id]=val;
+   }
+
+   /// Set the Velocity BC on a given ParGridFunction.
+   void SetEssVBC(mfem::ParGridFunction& pgf);
+
+   /// Extracts the true boundary doffs of the velocity
+   void SetEssTDofsV(mfem::Array<int>& ess_dofs);
+
+   /// Set the Velocity BC on a given true vector.
+   void SetEssTDofsV(mfem::Vector& v) const;
+
+   /// Forward solve with given RHS. x is the RHS vector.
+   /// The BC are set in the Mult operation.
+   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 in the MultTranspose operator.
+   void MultTranspose(const mfem::Vector &x, mfem::Vector &y) const override;
+
+   /// Return velocity grid function
+   mfem::ParGridFunction& GetVelocity()
+   {
+      vel.SetFromTrueDofs(sol.GetBlock(0));
+      return vel;
+   }
+
+   /// Return pressure grid function
+   mfem::ParGridFunction& GetPressure()
+   {
+      pre.SetFromTrueDofs(sol.GetBlock(1));
+      return pre;
+   }
+
+   /// Return velocity grid function
+   mfem::ParGridFunction& GetAdjVelocity()
+   {
+      avel.SetFromTrueDofs(adj.GetBlock(0));
+      return avel;
+   }
+
+   /// Return pressure grid function
+   mfem::ParGridFunction& GetAdjPressure()
+   {
+      apre.SetFromTrueDofs(adj.GetBlock(1));
+      return apre;
+   }
+
+   /// Return velocty FES
+   mfem::ParFiniteElementSpace* GetVelocitySpace() {return vfes;}
+
+   /// Return pressure FES
+   mfem::ParFiniteElementSpace* GetPressureSpace() {return pfes;}
+
+   /// Return the block offset for providing RHS vectors for
+   /// Mult and MultTranspose.
+   mfem::Array<int>& GetTrueBlockOffsets() {return block_true_offsets;}
+
+
+private:
+
+   int myrank;
+
+   bool zero_mean_press;
+
+   /// The parallel mesh.
+   mfem::ParMesh *pmesh = nullptr;
+
+   /// The order of the velocity  space.
+   int order;
+
+   /// linear system solvers parameters
+   real_t linear_atol;
+   real_t linear_rtol;
+   int  linear_iter;
+   int dim;
+
+   std::shared_ptr<mfem::Coefficient> visc; //viscosity
+   std::shared_ptr<mfem::Coefficient> brink; //Brinkman penalization
+
+   mfem::H1_FECollection* vfec; //velocity collections
+   mfem::FiniteElementCollection* pfec; //pressure collecation
+   mfem::ParFiniteElementSpace* vfes;
+   mfem::ParFiniteElementSpace* pfes;
+
+   std::unique_ptr<mfem::IterativeSolver> ls;
+   std::unique_ptr<mfem::IterativeSolver> prec;
+
+   //boundary conditions
+   std::map<int, std::shared_ptr<mfem::VectorCoefficient>> vel_bcs;
+
+   // holds the velocity constrained DOFs
+   mfem::Array<int> ess_tdofv;
+
+   // Volume force coefficient
+   std::shared_ptr<mfem::VectorCoefficient> vol_force;
+
+   mfem::Array<int> block_true_offsets;
+   int siz_u;
+   int siz_p;
+
+   mfem::ConstantCoefficient onecoeff;
+   mfem::ConstantCoefficient zerocoef;
+
+   mutable mfem::ParGridFunction vel; //velocity
+   mutable mfem::ParGridFunction pre; //pressure
+
+   mfem::ParGridFunction avel; //velocity
+   mfem::ParGridFunction apre; //pressure
+
+   std::unique_ptr<mfem::HypreParMatrix> A11;
+   std::unique_ptr<mfem::HypreParMatrix> A12;
+   std::unique_ptr<mfem::HypreParMatrix> A21;
+
+   std::unique_ptr<mfem::HypreParMatrix> A11e;
+   std::unique_ptr<mfem::HypreParMatrix> A12e;
+   std::unique_ptr<mfem::HypreParMatrix> A21e;
+
+   std::unique_ptr<mfem::ParBilinearForm> bf11;
+   std::unique_ptr<mfem::ParMixedBilinearForm> bf12;
+   std::unique_ptr<mfem::ParMixedBilinearForm> bf21;
+
+   std::unique_ptr<mfem::BlockOperator> bop;
+   mutable mfem::BlockVector rhs;
+   mutable mfem::BlockVector sol;
+   mutable mfem::BlockVector adj;
+
+   mfem::Vector V; //used for removing the mean pressure
+};
+
+class LSCStokesPrec:public mfem::IterativeSolver
+{
+public:
+   LSCStokesPrec(mfem::ParFiniteElementSpace* vfes_,
+                 mfem::ParFiniteElementSpace* pfes_,
+                 std::shared_ptr<mfem::Coefficient> visc_,
+                 std::shared_ptr<mfem::Coefficient> brink_,
+                 mfem::Array<int>& ess_vdofs_,
+                 const mfem::Operator* O11_,
+                 const mfem::Operator* O12_,
+                 const mfem::Operator* O21_,
+                 bool   zero_mean_press_=false)
+      :O11(O11_),O12(O12_),O21(O21_),zero_mean_press(zero_mean_press_)
+
+   {
+      // set the preconditioner for the upper block
+      mfem::ConstantCoefficient lambda(0.00);
+      std::unique_ptr<mfem::ParLORDiscretization>
+      lor_discr(new mfem::ParLORDiscretization(*vfes_));
+      mfem::ParFiniteElementSpace& vlor=lor_discr->GetParFESpace();
+      std::unique_ptr<mfem::ParBilinearForm>
+                             b11(new mfem::ParBilinearForm(&vlor));
+                             //b11(new mfem::ParBilinearForm(vfes_));
+      b11->AddDomainIntegrator(new mfem::ElasticityIntegrator(lambda,*visc_));
+      //b11->AddDomainIntegrator(new mfem::VectorDiffusionIntegrator(*visc_));
+      if (nullptr!=brink_.get())
+      {
+         b11->AddDomainIntegrator(new mfem::VectorMassIntegrator(*brink_));
+      }
+
+      b11->Assemble(0);
+      b11->Finalize(0);
+      A11.reset(b11->ParallelAssemble());
+      std::unique_ptr<mfem::HypreParMatrix> Ae(A11->EliminateRowsCols(ess_vdofs_));
+
+      std::cout<<"A11 assembled"<<std::endl;
+
+      amg11.reset(new mfem::HypreBoomerAMG());
+      if (mfem::Ordering::Type::byNODES==vfes_->GetOrdering())
+      {
+         int dim=vlor.GetParMesh()->Dimension();
+         amg11->SetSystemsOptions(dim,true);
+         //amg11->SetElasticityOptions(&vlor);
+      }
+      amg11->SetOperator(*A11);
+
+      cg11.reset(new mfem::CGSolver(vlor.GetComm()));
+      cg11->SetOperator(*O11);
+      cg11->SetPreconditioner(*amg11);
+      cg11->SetMaxIter(10);
+      cg11->SetRelTol(1e-12);
+      cg11->SetAbsTol(1e-12);
+      cg11->SetPrintLevel(0);
+
+      //assemble diagonal mass matrix on the velocity space
+      {
+         std::unique_ptr<mfem::ParBilinearForm>
+         q11(new mfem::ParBilinearForm(vfes_));
+         q11->AddDomainIntegrator(
+            new mfem::LumpedIntegrator(
+               //new mfem::VectorMassIntegrator(*brink_)));
+               new mfem::VectorMassIntegrator()));
+         q11->Assemble(0);
+         q11->Finalize(0);
+         Qv.reset(q11->ParallelAssemble());
+         amgv.reset(new mfem::HypreBoomerAMG());
+         amgv->SetOperator(*Qv);
+      }
+
+      const mfem::HypreParMatrix* m21=dynamic_cast<const mfem::HypreParMatrix*>(O21);
+      const mfem::HypreParMatrix* m12=dynamic_cast<const mfem::HypreParMatrix*>(O12);
+
+      if ((nullptr!=m12)&&(nullptr!=m21))
+      {
+
+         mfem::HypreParVector Sd(vfes_->GetComm(),
+                                 Qv->GetGlobalNumRows(),
+                                 Qv->GetRowStarts());
+         Qv->GetDiag(Sd);
+
+         mfem::HypreParMatrix T(*m12);
+         T.InvScaleRows(Sd);
+         A.reset(ParMult(m21,&T));
+      }
+      else
+      {
+         //Construct the discrete approximations for O12 and O21
+         std::unique_ptr<mfem::HypreParMatrix> A12, A21;
+         std::unique_ptr<mfem::ParMixedBilinearForm>
+         bf12(new mfem::ParMixedBilinearForm(pfes_, vfes_));
+         bf12->AddDomainIntegrator(
+            new mfem::TransposeIntegrator(
+               new mfem::VectorDivergenceIntegrator()));
+         bf12->Assemble(0);
+         bf12->Finalize(0);
+         A12.reset(bf12->ParallelAssemble());
+         A12->EliminateRows(ess_vdofs_);
+
+         std::unique_ptr<mfem::ParMixedBilinearForm>
+         bf21(new mfem::ParMixedBilinearForm(vfes_, pfes_));
+         bf21->AddDomainIntegrator(
+            new mfem::VectorDivergenceIntegrator());
+         bf21->Assemble(0);
+         bf21->Finalize(0);
+         A21.reset(bf21->ParallelAssemble());
+         std::unique_ptr<mfem::HypreParMatrix> A21e(A21->EliminateCols(ess_vdofs_));
+
+         mfem::HypreParVector Sd(vfes_->GetComm(),
+                                 Qv->GetGlobalNumRows(),
+                                 Qv->GetRowStarts());
+         Qv->GetDiag(Sd);
+
+         A12->InvScaleRows(Sd);
+         A.reset(ParMult(A21.get(),A12.get()));
+      }
+
+      amga.reset(new mfem::HypreBoomerAMG());
+      amga->SetOperator(*A);
+
+      if (zero_mean_press)
+      {
+         //use GMRES
+         cga.reset(new mfem::GMRESSolver(vfes_->GetComm()));
+         cga->SetOperator(*A);
+      }
+      else
+      {
+         //use CG
+         cga.reset(new mfem::CGSolver(vfes_->GetComm()));
+         cga->SetOperator(*A);
+      }
+
+      cga->SetPreconditioner(*amga);
+      cga->SetPrintLevel(0);
+      cga->SetMaxIter(20);
+      cga->SetRelTol(1e-12);
+      cga->SetAbsTol(1e-12);
+
+      siz_u=O11->NumRows();
+      siz_p=O21->NumRows();
+
+      block_true_offsets.SetSize(3);
+      block_true_offsets[0] = 0;
+      block_true_offsets[1] = siz_u;
+      block_true_offsets[2] = siz_p;
+      block_true_offsets.PartialSum();
+      //set the width and the height of the operator
+      this->width=  block_true_offsets[2];
+      this->height= block_true_offsets[2];
+
+      v1.SetSize(siz_p); v1=0.0;
+      v2.SetSize(siz_u); v2=0.0;
+      v3.SetSize(siz_u); v3=0.0;
+      v4.SetSize(siz_p); v4=0.0;
+
+      myrank=vfes_->GetMyRank();
+
+
+   }
+
+   /// Operator application
+   virtual
+   void Mult (const mfem::Vector & x, mfem::Vector & y) const override
+   {
+      mfem::BlockVector xb,yb;
+
+      cga->SetMaxIter(mfem::IterativeSolver::max_iter);
+      cga->SetAbsTol(IterativeSolver::abs_tol);
+      cga->SetRelTol(IterativeSolver::rel_tol);
+      //cga->iterative_mode=true;
+
+      cg11->SetMaxIter(mfem::IterativeSolver::max_iter);
+      cg11->SetAbsTol(IterativeSolver::abs_tol);
+      cg11->SetRelTol(IterativeSolver::rel_tol);
+      //cg11->iterative_mode=true;
+
+
+      xb.Update(const_cast<mfem::Vector&>(x), block_true_offsets);
+      yb.Update(y, block_true_offsets);
+
+      if(0==myrank){std::cout<<"Schur complement solve";}
+      cga->Mult(xb.GetBlock(1),v1);
+      O12->Mult(v1,v2);
+      amgv->Mult(v2,v3);
+      O11->Mult(v3,v2);
+      amgv->Mult(v2,v3);
+      O21->Mult(v3,v4);
+      cga->Mult(v4,yb.GetBlock(1));
+      yb.GetBlock(1).Neg();
+
+      //construct modification of the rhs for block 0
+      O12->Mult(yb.GetBlock(1),v2);
+      add(xb.GetBlock(0), -1, v2, v3);
+
+      //multiply the upper block
+      if(0==myrank){std::cout<<"Upper block solve";}
+      cg11->Mult(v3,yb.GetBlock(0));
+      //amg11->Mult(v3,yb.GetBlock(0));
+   }
+
+private:
+
+   mutable mfem::Vector v1;
+   mutable mfem::Vector v2;
+   mutable mfem::Vector v3;
+   mutable mfem::Vector v4;
+
+   const mfem::Operator* O11;
+   const mfem::Operator* O12;
+   const mfem::Operator* O21;
+
+   std::unique_ptr<mfem::HypreParMatrix> A11;
+   std::unique_ptr<mfem::HypreBoomerAMG> amg11;
+   std::unique_ptr<mfem::CGSolver> cg11;
+
+   std::unique_ptr<mfem::HypreParMatrix> A;
+   std::unique_ptr<mfem::HypreBoomerAMG> amga;
+   std::unique_ptr<mfem::IterativeSolver> cga;
+
+
+   int siz_u;
+   int siz_p;
+
+   mfem::Array<int> block_true_offsets;
+
+   std::unique_ptr<mfem::HypreParMatrix> Qv;
+   std::unique_ptr<mfem::HypreBoomerAMG> amgv;
+
+   bool zero_mean_press;
+   int myrank;
+};

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/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/test_grain.cpp

@@ -0,0 +1,156 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+#include "mtop_solvers.hpp"
+#include "grain_reader.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 = "./mini_flow2d_ball.msh";
+   int order = 1;
+   bool static_cond = false;
+   int ser_ref_levels = 0;
+   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);
+
+   double meshOffsetX = -12.0;
+   double meshOffsetY = -12.0;
+   double meshOffsetZ = -0.5;
+
+
+   double Lx = 24.0; double Ly = 24.0; double Lz = 150.5;
+   int NX = 48;      int NY = 48;      int NZ = 301;
+   mfem::Mesh mesh = mfem::Mesh::MakeCartesian3D(NX, NY, NZ, mfem::Element::HEXAHEDRON, Lx, Ly, Lz, true);
+   int dim = mesh.Dimension();
+
+   int tNumVertices  = mesh.GetNV();
+   for (int i = 0; i < tNumVertices; ++i)
+   {
+      double * Coords = mesh.GetVertex(i);
+
+      Coords[ 0 ] = Coords[ 0 ] + meshOffsetX;
+      Coords[ 1 ] = Coords[ 1 ] + meshOffsetY;
+      Coords[ 2 ] = Coords[ 2 ] + meshOffsetZ;
+   }
+
+
+
+   // 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.
+   {
+      for (int l = 0; l < ser_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();
+      }
+   }
+
+   std::cout<<"My rank="<<pmesh.GetMyRank()<<std::endl;
+   std::string tStringWeight = "./singleparticlesdiameter.txt";
+   ::mfem::H1_FECollection FECol_H1(order, dim);
+   ::mfem::ParFiniteElementSpace FESpace_H1(&pmesh, &FECol_H1, 1, mfem::Ordering::byNODES);
+
+   GrainReader grain( &pmesh, tStringWeight );
+   grain.computeGridFunction( FESpace_H1);
+   ::mfem::ParGridFunction grainLSField =  grain.getGrainGridFunction();
+
+   //dump the solution
+   {
+      ParaViewDataCollection paraview_dc("grain_test", &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("level_set",&grainLSField);
+
+      paraview_dc.Save();
+   }
+
+
+   MPI::Finalize();
+   return 0;
+}
+
+

miniapps/mtop/test_stokes.cpp

@@ -0,0 +1,203 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+#include "mtop_solvers.hpp"
+
+using namespace mfem;
+
+class BrinkCoeff :public Coefficient
+{
+public:
+    BrinkCoeff(real_t penal_=10.0):penalty(penal_)
+    {
+    }
+
+    virtual real_t Eval(ElementTransformation &T,
+                        const IntegrationPoint &ip) override
+    {
+        real_t x[3];
+        Vector transip(x, 3);
+        T.Transform(ip, transip);
+
+        real_t c[3]={0.2,0.2,0.0};
+        real_t r=0.0;
+        real_t d;
+        for(int i=0;i<3;i++){
+            d=x[i]-c[i];
+            r=r+d*d;
+        }
+
+        r=sqrt(r);
+        if(r>0.05){return 0.0;}
+        else{ return penalty;}
+    }
+
+private:
+    real_t penalty;
+};
+
+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 = "./dfg_bench_flow_tri.msh";
+   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(1000./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();
+      }
+   }
+
+   std::cout<<"My rank="<<pmesh.GetMyRank()<<std::endl;
+
+   StokesSolver* solver=new StokesSolver(&pmesh,2);
+
+
+   mfem::Vector bci(dim); bci=1.0; bci(1)=0.0;
+   mfem::Vector zvi(dim); zvi=0.0;
+   std::shared_ptr<VectorCoefficient> cvci;
+   cvci.reset(new VectorConstantCoefficient(bci));
+   std::shared_ptr<VectorCoefficient> zvci;
+   zvci.reset(new VectorConstantCoefficient(zvi));
+
+   solver->AddVelocityBC(1,cvci);
+   solver->AddVelocityBC(2,cvci);
+   solver->AddVelocityBC(3,cvci);
+   //solver->AddVelocityBC(4,cvci);
+   //solver->AddVelocityBC(5,zvci);
+
+   std::shared_ptr<Coefficient> brink;
+   brink.reset(new BrinkCoeff(1000.0));
+
+   solver->SetBrink(brink);
+
+   ParGridFunction pg(solver->GetVelocitySpace()); pg=0.0;
+   ParGridFunction ng(solver->GetVelocitySpace()); ng=0.0;
+   ng.SetTrueVector();
+   solver->SetEssVBC(pg);
+   Vector pgv(ng.GetTrueVector());
+
+   solver->SetEssTDofsV(pgv);
+
+   ng.SetFromTrueDofs(pgv);
+
+   //solver->SetZeroMeanPressure(true);
+   solver->SetLinearSolver(1e-8,1e-12,550);
+
+   solver->Assemble();
+   solver->FSolve();
+
+   //dump the solution
+   {
+      ParGridFunction& vel=solver->GetVelocity();
+      ParGridFunction& pre=solver->GetPressure();
+
+      ParaViewDataCollection paraview_dc("stokes_flow", &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("vel",&vel);
+      paraview_dc.RegisterField("pres",&pre);
+      paraview_dc.RegisterField("pg",&pg);
+      paraview_dc.RegisterField("ng",&ng);
+      paraview_dc.Save();
+   }
+
+   delete solver;
+
+   MPI::Finalize();
+   return 0;
+}
+
+