Add unit test for operator ownership in sum and product operators

SumOperator, ProductOperator, and TripleProductOperator use Handle<T> instead of
raw pointers. Explicit ownership flags and destructors can be removed. The
classes now have proper copy and move semantics (rule of zero).

Retain the constructors with explicit ownership flags for backwards
compatibility.

.gitignore

@@ -8,6 +8,7 @@
 
 # Object and library files
 *.o
+*.o.tmp
 /libmfem.*
 /miniapps/common/libmfem-common.*
 

CHANGELOG

@@ -75,6 +75,11 @@ Miscellaneous
 
 - Added support for custom interpolation procedure in FindPointsGSLIB.
 
+- `FiniteElementSpace` has new methods to directly set prolongation and
+  restriction operators to arbitrary sparse matrices.
+
+- There are new convenience constructors for NURBS patches and knot vectors.
+
 API changes
 -----------
 - API change: in class GridFunction, 'fec' was renamed to 'fec_owned'.

fem/bilininteg.cpp

@@ -855,6 +855,34 @@ const IntegrationRule &GradientIntegrator::GetRule(const FiniteElement
 }
 
 
+DiffusionIntegrator::DiffusionIntegrator(const IntegrationRule *ir)
+   : BilinearFormIntegrator(ir),
+     Q(nullptr), VQ(nullptr), MQ(nullptr), maps(nullptr), geom(nullptr)
+{
+   static Kernels kernels;
+}
+
+DiffusionIntegrator::DiffusionIntegrator(Coefficient &q,
+                                         const IntegrationRule *ir)
+   : DiffusionIntegrator(ir)
+{
+   Q = &q;
+}
+
+DiffusionIntegrator::DiffusionIntegrator(VectorCoefficient &q,
+                                         const IntegrationRule *ir)
+   : DiffusionIntegrator(ir)
+{
+   VQ = &q;
+}
+
+DiffusionIntegrator::DiffusionIntegrator(MatrixCoefficient &q,
+                                         const IntegrationRule *ir)
+   : DiffusionIntegrator(ir)
+{
+   MQ = &q;
+}
+
 void DiffusionIntegrator::AssembleElementMatrix
 ( const FiniteElement &el, ElementTransformation &Trans,
   DenseMatrix &elmat )
@@ -1310,6 +1338,17 @@ const IntegrationRule &DiffusionIntegrator::GetRule(
    return IntRules.Get(trial_fe.GetGeomType(), order);
 }
 
+MassIntegrator::MassIntegrator(const IntegrationRule *ir)
+   : BilinearFormIntegrator(ir), Q(nullptr), maps(nullptr), geom(nullptr)
+{
+   static Kernels kernels;
+}
+
+MassIntegrator::MassIntegrator(Coefficient &q, const IntegrationRule *ir)
+   : MassIntegrator(ir)
+{
+   Q = &q;
+}
 
 void MassIntegrator::AssembleElementMatrix
 ( const FiniteElement &el, ElementTransformation &Trans,

fem/bilininteg.hpp

@@ -2156,7 +2156,7 @@ public:
 
    MFEM_REGISTER_KERNELS(ApplyPAKernels, ApplyKernelType, (int, int, int));
    MFEM_REGISTER_KERNELS(DiagonalPAKernels, DiagonalKernelType, (int, int, int));
-   static struct Kernels { Kernels(); } kernels;
+   struct Kernels { Kernels(); };
 
 protected:
    Coefficient *Q;
@@ -2234,26 +2234,16 @@ private:
 
 public:
    /// Construct a diffusion integrator with coefficient Q = 1
-   DiffusionIntegrator(const IntegrationRule *ir = nullptr)
-      : BilinearFormIntegrator(ir),
-        Q(NULL), VQ(NULL), MQ(NULL), maps(NULL), geom(NULL) { }
+   DiffusionIntegrator(const IntegrationRule *ir = nullptr);
 
    /// Construct a diffusion integrator with a scalar coefficient q
-   DiffusionIntegrator(Coefficient &q, const IntegrationRule *ir = nullptr)
-      : BilinearFormIntegrator(ir),
-        Q(&q), VQ(NULL), MQ(NULL), maps(NULL), geom(NULL) { }
+   DiffusionIntegrator(Coefficient &q, const IntegrationRule *ir = nullptr);
 
    /// Construct a diffusion integrator with a vector coefficient q
-   DiffusionIntegrator(VectorCoefficient &q,
-                       const IntegrationRule *ir = nullptr)
-      : BilinearFormIntegrator(ir),
-        Q(NULL), VQ(&q), MQ(NULL), maps(NULL), geom(NULL) { }
+   DiffusionIntegrator(VectorCoefficient &q, const IntegrationRule *ir = nullptr);
 
    /// Construct a diffusion integrator with a matrix coefficient q
-   DiffusionIntegrator(MatrixCoefficient &q,
-                       const IntegrationRule *ir = nullptr)
-      : BilinearFormIntegrator(ir),
-        Q(NULL), VQ(NULL), MQ(&q), maps(NULL), geom(NULL) { }
+   DiffusionIntegrator(MatrixCoefficient &q, const IntegrationRule *ir = nullptr);
 
    /** Given a particular Finite Element computes the element stiffness matrix
        elmat. */
@@ -2356,15 +2346,13 @@ public:
 
    MFEM_REGISTER_KERNELS(ApplyPAKernels, ApplyKernelType, (int, int, int));
    MFEM_REGISTER_KERNELS(DiagonalPAKernels, DiagonalKernelType, (int, int, int));
-   static struct Kernels { Kernels(); } kernels;
+   struct Kernels { Kernels(); };
 
 public:
-   MassIntegrator(const IntegrationRule *ir = NULL)
-      : BilinearFormIntegrator(ir), Q(NULL), maps(NULL), geom(NULL) { }
+   MassIntegrator(const IntegrationRule *ir = nullptr);
 
    /// Construct a mass integrator with coefficient q
-   MassIntegrator(Coefficient &q, const IntegrationRule *ir = NULL)
-      : BilinearFormIntegrator(ir), Q(&q), maps(NULL), geom(NULL) { }
+   MassIntegrator(Coefficient &q, const IntegrationRule *ir = NULL);
 
    /** Given a particular Finite Element computes the element mass matrix
        elmat. */

fem/conduitdatacollection.cpp

@@ -997,8 +997,7 @@ std::string
 ConduitDataCollection::MeshFilePattern(const std::string &relay_protocol)
 {
    std::ostringstream oss;
-   oss << prefix_path
-       << name
+   oss << name
        << "_"
        << to_padded_string(cycle, pad_digits_cycle)
        << "/domain_%0"

fem/fe/face_map_utils.cpp

@@ -41,8 +41,8 @@ void FillFaceMap(const int n_face_dofs_per_component,
                  const std::vector<int> &n_dofs_per_dim,
                  Array<int> &face_map)
 {
-   const int n_components = offsets.size();
-   const int face_dim = strides.size() / n_components;
+   const int n_components = static_cast<int>(offsets.size());
+   const int face_dim = static_cast<int>(strides.size()) / n_components;
    for (int comp = 0; comp < n_components; ++comp)
    {
       const int offset = offsets[comp];

fem/fespace.cpp

@@ -146,6 +146,43 @@ void FiniteElementSpace::CopyProlongationAndRestriction(
    delete perm_mat_tr;
 }
 
+void FiniteElementSpace::SetProlongation(const SparseMatrix& p)
+{
+#ifdef MFEM_USE_MPI
+   MFEM_VERIFY(dynamic_cast<const ParFiniteElementSpace*>(this) == NULL,
+               "Attempting to set serial prolongation operator for "
+               "parallel finite element space.");
+#endif
+
+   if (!cP)
+   {
+      cP = std::unique_ptr<SparseMatrix>(new SparseMatrix(p));
+   }
+   else
+   {
+      *cP = p;
+   }
+   cP_is_set = true;
+}
+
+void FiniteElementSpace::SetRestriction(const SparseMatrix& r)
+{
+#ifdef MFEM_USE_MPI
+   MFEM_VERIFY(dynamic_cast<const ParFiniteElementSpace*>(this) == NULL,
+               "Attempting to set serial restriction operator for "
+               "parallel finite element space.");
+#endif
+
+   if (!cR)
+   {
+      cR = std::unique_ptr<SparseMatrix>(new SparseMatrix(r));
+   }
+   else
+   {
+      *cR = r;
+   }
+}
+
 void FiniteElementSpace::SetElementOrder(int i, int p)
 {
    MFEM_VERIFY(mesh_sequence == mesh->GetSequence(),

fem/fespace.hpp

@@ -587,6 +587,14 @@ public:
    bool Conforming() const { return mesh->Conforming() && cP == NULL; }
    bool Nonconforming() const { return mesh->Nonconforming() || cP != NULL; }
 
+   /** Set the prolongation operator of the space to an arbitrary sparse matrix,
+       creating a copy of the argument. */
+   void SetProlongation(const SparseMatrix& p);
+
+   /** Set the restriction operator of the space to an arbitrary sparse matrix,
+       creating a copy of the argument. */
+   void SetRestriction(const SparseMatrix& r);
+
    /// Sets the order of the i'th finite element.
    /** By default, all elements are assumed to be of fec->GetOrder(). Once
        SetElementOrder is called, the space becomes a variable order space. */

fem/gslib.cpp

@@ -55,7 +55,7 @@ FindPointsGSLIB::FindPointsGSLIB()
    gsl_comm = new gslib::comm;
    cr       = new gslib::crystal;
 #ifdef MFEM_USE_MPI
-   int initialized;
+   int initialized = 0;
    MPI_Initialized(&initialized);
    if (!initialized) { MPI_Init(NULL, NULL); }
    MPI_Comm comm = MPI_COMM_WORLD;

fem/integ/bilininteg_diffusion_kernels.cpp

@@ -16,7 +16,6 @@ namespace mfem
 
 // PA Diffusion Integrator
 
-DiffusionIntegrator::Kernels DiffusionIntegrator::kernels;
 DiffusionIntegrator::Kernels::Kernels()
 {
    // 2D

fem/integ/bilininteg_diffusion_patch.cpp

@@ -1039,7 +1039,7 @@ void DiffusionIntegrator::AssemblePatchMatrix_reducedQuadrature(
       for (int zquad = 0; zquad<2; ++zquad)
       {
          // Reduced quadrature in z
-         const int nwz = rid(zquad,2,patch)[jdz].size();
+         const int nwz = static_cast<int>(rid(zquad,2,patch)[jdz].size());
          for (int irz=0; irz < nwz; ++irz)
          {
             const int qz = rid(zquad,2,patch)[jdz][irz] + minD[2][jdz];
@@ -1062,7 +1062,7 @@ void DiffusionIntegrator::AssemblePatchMatrix_reducedQuadrature(
             for (int yquad = 0; yquad<2; ++yquad)
             {
                // Reduced quadrature in y
-               const int nwy = rid(yquad,1,patch)[jdy].size();
+               const int nwy = static_cast<int>(rid(yquad,1,patch)[jdy].size());
                for (int iry=0; iry < nwy; ++iry)
                {
                   const int qy = rid(yquad,1,patch)[jdy][iry] + minD[1][jdy];
@@ -1082,7 +1082,7 @@ void DiffusionIntegrator::AssemblePatchMatrix_reducedQuadrature(
                   // Reduced quadrature in x
                   for (int xquad=0; xquad<2; ++xquad)
                   {
-                     const int nwx = rid(xquad,0,patch)[jdx].size();
+                     const int nwx = static_cast<int>(rid(xquad,0,patch)[jdx].size());
                      for (int irx=0; irx < nwx; ++irx)
                      {
                         const int qx = rid(xquad,0,patch)[jdx][irx] + minD[0][jdx];
@@ -1117,7 +1117,7 @@ void DiffusionIntegrator::AssemblePatchMatrix_reducedQuadrature(
                   }
 
                   // 00 terms
-                  const int nw = rid(0,0,patch)[jdx].size();
+                  const int nw = static_cast<int>(rid(0,0,patch)[jdx].size());
                   for (int irx=0; irx < nw; ++irx)
                   {
                      const int qx = rid(0,0,patch)[jdx][irx] + minD[0][jdx];
@@ -1140,7 +1140,7 @@ void DiffusionIntegrator::AssemblePatchMatrix_reducedQuadrature(
                   }
 
                   // 11 terms
-                  const int nw11 = rid(1,0,patch)[jdx].size();
+                  const int nw11 = static_cast<int>(rid(1,0,patch)[jdx].size());
 
                   for (int irx=0; irx < nw11; ++irx)
                   {

fem/integ/bilininteg_mass_kernels.cpp

@@ -14,7 +14,6 @@
 namespace mfem
 {
 
-MassIntegrator::Kernels MassIntegrator::kernels;
 MassIntegrator::Kernels::Kernels()
 {
    // 2D

fem/intrules.cpp

@@ -1908,7 +1908,7 @@ IntegrationRule& NURBSMeshRules::GetElementRule(const int elem,
          }
       }
 
-      npd[d] = el[d].size() / 2;
+      npd[d] = static_cast<int>(el[d].size() / 2);
       np *= npd[d];
    }
 

fem/kernel_dispatch.hpp

@@ -78,9 +78,9 @@ namespace mfem
       const char *kernel_name = MFEM_KERNEL_NAME(KernelName);                  \
       using KernelSignature = KernelType;                                      \
       template <MFEM_PARAM_LIST P3>                                            \
-      static KernelSignature Kernel();                                         \
-      static KernelSignature Fallback(MFEM_PARAM_LIST P1);                     \
-      static KernelName &Get()                                                 \
+      static MFEM_EXPORT KernelSignature Kernel();                             \
+      static MFEM_EXPORT KernelSignature Fallback(MFEM_PARAM_LIST P1);         \
+      static MFEM_EXPORT KernelName &Get()                                     \
       { static KernelName table; return table;}                                \
    }
 
@@ -126,9 +126,9 @@ class KernelDispatchTable<Kernels,
          internal::KernelTypeList<Params...>,
          internal::KernelTypeList<OptParams...>>
 {
-   std::unordered_map<std::tuple<Params...>,
-       Signature,
-       KernelDispatchKeyHash<Params...>> table;
+   using TableType = std::unordered_map<std::tuple<Params...>,
+         Signature, KernelDispatchKeyHash<Params...>>;
+   TableType table;
 
 public:
    /// @brief Run the kernel with the given dispatch parameters and arguments.
@@ -176,6 +176,12 @@ public:
          }
       };
    };
+
+   /// Return the dispatch map table
+   static const TableType &GetDispatchTable()
+   {
+      return Kernels::Get().table;
+   }
 };
 
 }

fem/lininteg.hpp

@@ -676,7 +676,7 @@ public:
       int myid;
       MPI_Comm_rank(comm, &myid);
 
-      int seed = (seed_ > 0) ? seed_ + myid : time(0) + myid;
+      int seed = (seed_ > 0) ? seed_ + myid : (int)time(0) + myid;
       SetSeed(seed);
    }
 #else

fem/multigrid.hpp

@@ -16,7 +16,7 @@
 #include "bilinearform.hpp"
 
 #include "../linalg/operator.hpp"
-#include "../linalg/handle.hpp"
+#include "../linalg/op_handle.hpp"
 
 namespace mfem
 {

fem/pfespace.cpp

@@ -1984,7 +1984,7 @@ struct PMatrixRow
 
    void write(std::ostream &os, real_t sign) const
    {
-      bin_io::write<int>(os, elems.size());
+      bin_io::write<int>(os, static_cast<int>(elems.size()));
       for (unsigned i = 0; i < elems.size(); i++)
       {
          const PMatrixElement &e = elems[i];
@@ -2074,7 +2074,7 @@ void NeighborRowMessage::Encode(int rank)
    }
 
    Array<GroupId> all_group_ids;
-   all_group_ids.Reserve(rows.size());
+   all_group_ids.Reserve(static_cast<int>(rows.size()));
    for (int i = 0; i < 3; i++)
    {
       all_group_ids.Append(group_ids[i]);
@@ -2833,7 +2833,7 @@ HypreParMatrix* ParFiniteElementSpace
    }
 
    // create offd column mapping
-   HYPRE_BigInt *cmap = Memory<HYPRE_BigInt>(col_map.size());
+   HYPRE_BigInt *cmap = Memory<HYPRE_BigInt>(static_cast<int>(col_map.size()));
    int offd_col = 0;
    for (auto it = col_map.begin(); it != col_map.end(); ++it)
    {
@@ -2893,7 +2893,7 @@ HypreParMatrix* ParFiniteElementSpace
                              row_starts.GetData(), col_starts.GetData(),
                              I_diag, J_diag, A_diag,
                              I_offd, J_offd, A_offd,
-                             col_map.size(), cmap);
+                             static_cast<HYPRE_Int>(col_map.size()), cmap);
 }
 
 template <typename int_type>
@@ -3119,7 +3119,7 @@ ParFiniteElementSpace::ParallelDerefinementMatrix(int old_ndofs,
             msg.dofs[i] = old_offset + dofs[i];
          }
 
-         MPI_Isend(&msg.dofs[0], msg.dofs.size(), HYPRE_MPI_BIG_INT,
+         MPI_Isend(&msg.dofs[0], static_cast<int>(msg.dofs.size()), HYPRE_MPI_BIG_INT,
                    coarse_rank, 291, MyComm, &msg.request);
       }
       else if (coarse_rank == MyRank && fine_rank != MyRank)
@@ -3240,7 +3240,7 @@ ParFiniteElementSpace::ParallelDerefinementMatrix(int old_ndofs,
                   {
                      if (row[j] == 0.0) { continue; } // NOTE: lR thresholded
                      int &lcol = col_map[remote_dofs[j]];
-                     if (!lcol) { lcol = col_map.size(); }
+                     if (!lcol) { lcol = static_cast<int>(col_map.size()); }
                      offd->_Set_(m, lcol-1, row[j]);
                   }
                   mark[m] = 1;
@@ -3252,7 +3252,7 @@ ParFiniteElementSpace::ParallelDerefinementMatrix(int old_ndofs,
 
    messages.clear();
    offd->Finalize(0);
-   offd->SetWidth(col_map.size());
+   offd->SetWidth(static_cast<int>(col_map.size()));
 
    // create offd column mapping for use by hypre
    HYPRE_BigInt *cmap = Memory<HYPRE_BigInt>(offd->Width());

fem/pgridfunc.cpp

@@ -982,70 +982,65 @@ void ParGridFunction::SaveAsSerial(const char *fname, int precision,
    MPI_Barrier(pmesh->GetComm());
 }
 
-GridFunction ParGridFunction::GetSerialGridFunction(int save_rank,
-                                                    Mesh &serial_mesh) const
+GridFunction ParGridFunction::GetSerialGridFunction(
+   int save_rank, FiniteElementSpace &serial_fes) const
 {
    ParFiniteElementSpace *pfespace = ParFESpace();
    ParMesh *pmesh = pfespace->GetParMesh();
 
-   int vdim = pfespace->GetVDim();
-   auto *fec_serial = FiniteElementCollection::New(pfespace->FEColl()->Name());
-   auto *fespace_serial = new FiniteElementSpace(&serial_mesh,
-                                                 fec_serial,
-                                                 vdim,
-                                                 pfespace->GetOrdering());
+   GridFunction serial_gf(&serial_fes);
 
-   GridFunction gf_serial(fespace_serial);
-   gf_serial.MakeOwner(fec_serial);
    Array<real_t> vals;
    Array<int> dofs;
    MPI_Status status;
-   int n_send_recv;
 
-   int my_rank = pmesh->GetMyRank(),
-       nranks = pmesh->GetNRanks();
-   MPI_Comm my_comm = pmesh->GetComm();
+   const int vdim = pfespace->GetVDim();
 
-   int elem_count = 0; // To keep track of element count in serial mesh
+   const int my_rank = pmesh->GetMyRank();
+   const int nranks = pmesh->GetNRanks();
+   MPI_Comm comm = pmesh->GetComm();
 
    if (my_rank == save_rank)
    {
+      int elem_count = 0; // To keep track of element count in serial mesh
+
       Vector nodeval;
       for (int e = 0; e < pmesh->GetNE(); e++)
       {
          GetElementDofValues(e, nodeval);
-         fespace_serial->GetElementVDofs(elem_count++, dofs);
-         gf_serial.SetSubVector(dofs, nodeval);
+         serial_fes.GetElementVDofs(elem_count++, dofs);
+         serial_gf.SetSubVector(dofs, nodeval);
       }
 
       for (int p = 0; p < nranks; p++)
       {
          if (p == save_rank) { continue; }
-         MPI_Recv(&n_send_recv, 1, MPI_INT, p, 448, my_comm, &status);
+         int n_send_recv;
+         MPI_Recv(&n_send_recv, 1, MPI_INT, p, 448, comm, &status);
          vals.SetSize(n_send_recv);
          if (n_send_recv)
          {
-            MPI_Recv(&vals[0], n_send_recv, MPITypeMap<real_t>::mpi_type, p, 449, my_comm,
+            MPI_Recv(&vals[0], n_send_recv, MPITypeMap<real_t>::mpi_type, p, 449, comm,
                      &status);
          }
          for (int i = 0; i < n_send_recv; )
          {
-            fespace_serial->GetElementVDofs(elem_count++, dofs);
-            gf_serial.SetSubVector(dofs, &vals[i]);
+            serial_fes.GetElementVDofs(elem_count++, dofs);
+            serial_gf.SetSubVector(dofs, &vals[i]);
             i += dofs.Size();
          }
       }
    } // my_rank == save_rank
    else
    {
-      n_send_recv = 0;
+      int n_send_recv = 0;
       Vector nodeval;
       for (int e = 0; e < pmesh->GetNE(); e++)
       {
          const FiniteElement *fe = pfespace->GetFE(e);
          n_send_recv += vdim*fe->GetDof();
       }
-      MPI_Send(&n_send_recv, 1, MPI_INT, save_rank, 448, my_comm);
+      MPI_Send(&n_send_recv, 1, MPI_INT, save_rank, 448, comm);
       vals.Reserve(n_send_recv);
       vals.SetSize(0);
       for (int e = 0; e < pmesh->GetNE(); e++)
@@ -1059,12 +1054,24 @@ GridFunction ParGridFunction::GetSerialGridFunction(int save_rank,
       if (n_send_recv)
       {
          MPI_Send(&vals[0], n_send_recv, MPITypeMap<real_t>::mpi_type, save_rank, 449,
-                  my_comm);
+                  comm);
       }
    }
 
-   MPI_Barrier(my_comm);
-   return gf_serial;
+   return serial_gf;
+}
+
+GridFunction ParGridFunction::GetSerialGridFunction(int save_rank,
+                                                    Mesh &serial_mesh) const
+{
+   auto *serial_fec = pfes->FEColl()->Clone(pfes->FEColl()->GetOrder());
+   auto *serial_fes = new FiniteElementSpace(&serial_mesh,
+                                             serial_fec,
+                                             pfes->GetVDim(),
+                                             pfes->GetOrdering());
+   GridFunction serial_gf = GetSerialGridFunction(save_rank, *serial_fes);
+   serial_gf.MakeOwner(serial_fec); // Also assumes ownership of serial_fes
+   return serial_gf;
 }
 
 #ifdef MFEM_USE_ADIOS2

fem/pgridfunc.hpp

@@ -434,13 +434,29 @@ public:
    /// be used for ASCII output.
    void Save(const char *fname, int precision=16) const override;
 
-   /// Returns a GridFunction on MPI rank @a save_rank that does not have any
-   /// duplication of vertices/nodes at processor boundaries.
-   /// serial_mesh is obtained using ParMesh::GetSerialMesh(save_rank).
-   /// Note that the @ save_rank argument must match for the
-   /// ParMesh::GetSerialMesh and GetSerialGridFunction method.
+   /// @brief Returns a GridFunction on MPI rank @a save_rank that does not have
+   /// any duplication of vertices/nodes at processor boundaries.
+   ///
+   /// The @a serial_mesh is obtained using ParMesh::GetSerialMesh. Note that
+   /// the @a save_rank must be the same as that used in ParMesh::GetSerialMesh.
+   ///
+   /// @note The returned GridFunction will own the newly created
+   /// FiniteElementCollection and FiniteElementSpace objects.
    GridFunction GetSerialGridFunction(int save_rank, Mesh &serial_mesh) const;
 
+   /// @brief Returns a GridFunction on MPI rank @a save_rank that does not have
+   /// any duplication of vertices/nodes at processor boundaries.
+   ///
+   /// The given @a serial_fes must be defined on the mesh returned by
+   /// ParMesh::GetSerialMesh (with @a save_rank ranks), for example using the
+   /// space belonging to the GridFunction obtained from @ref
+   /// ParGridFunction::GetSerialGridFunction(int,Mesh &) const.
+   ///
+   /// @note The returned GridFunction does not assume ownership of @a
+   /// serial_fes.
+   GridFunction GetSerialGridFunction(
+      int save_rank, FiniteElementSpace &serial_fes) const;
+
    /// Write the serial GridFunction a single file (written using MPI rank 0).
    /// The given @a precision will be used for ASCII output.
    void SaveAsSerial(const char *fname, int precision=16, int save_rank=0) const;

fem/quadinterpolator.cpp

@@ -30,26 +30,27 @@ void InitEvalKernels();
 void InitDetKernels();
 template <bool P> void InitGradByNodesKernels();
 template <bool P> void InitGradByVDimKernels();
-}
-}
-
-QuadratureInterpolator::Kernels QuadratureInterpolator::kernels;
-QuadratureInterpolator::Kernels::Kernels()
+struct Kernels
 {
-   using namespace internal::quadrature_interpolator;
+   Kernels()
+   {
+      using namespace internal::quadrature_interpolator;
 
-   InitEvalByNodesKernels();
-   InitEvalByVDimKernels();
-   // Non-phys grad kernels
-   InitGradByNodesKernels<false>();
-   InitGradByVDimKernels<false>();
-   // Phys grad kernels
-   InitGradByNodesKernels<true>();
-   InitGradByVDimKernels<true>();
-   // Determinants
-   InitDetKernels();
-   // Non-tensor
-   InitEvalKernels();
+      InitEvalByNodesKernels();
+      InitEvalByVDimKernels();
+      // Non-phys grad kernels
+      InitGradByNodesKernels<false>();
+      InitGradByVDimKernels<false>();
+      // Phys grad kernels
+      InitGradByNodesKernels<true>();
+      InitGradByVDimKernels<true>();
+      // Determinants
+      InitDetKernels();
+      // Non-tensor
+      InitEvalKernels();
+   }
+};
+}
 }
 
 QuadratureInterpolator::QuadratureInterpolator(const FiniteElementSpace &fes,
@@ -61,6 +62,8 @@ QuadratureInterpolator::QuadratureInterpolator(const FiniteElementSpace &fes,
    q_layout(QVectorLayout::byNODES),
    use_tensor_products(UsesTensorBasis(fes))
 {
+   static internal::quadrature_interpolator::Kernels kernels;
+
    d_buffer.UseDevice(true);
    if (fespace->GetNE() == 0) { return; }
    const FiniteElement *fe = fespace->GetFE(0);

fem/quadinterpolator.hpp

@@ -158,8 +158,6 @@ public:
    MFEM_REGISTER_KERNELS(EvalKernels, EvalKernelType, (int, int, int, int));
    MFEM_REGISTER_KERNELS(CollocatedGradKernels, CollocatedGradKernelType,
                          (int, QVectorLayout, bool, int, int), (int));
-
-   static struct Kernels { Kernels(); } kernels;
 };
 
 }

fem/restriction.cpp

@@ -1631,7 +1631,7 @@ void InterpolationManager::LinearizeInterpolatorMapIntoVector()
    const FiniteElement *trace_fe =
       fes.GetTraceElement(0, fes.GetMesh()->GetFaceGeometry(0));
    const int face_dofs = trace_fe->GetDof();
-   const int nc_size = interp_map.size();
+   const int nc_size = static_cast<int>(interp_map.size());
    MFEM_VERIFY(nc_cpt==nc_size, "Unexpected number of interpolators.");
    interpolators.SetSize(face_dofs*face_dofs*nc_size);
    auto d_interp = Reshape(interpolators.HostWrite(),face_dofs,face_dofs,nc_size);

fem/tmop.cpp

@@ -1926,7 +1926,6 @@ void DiscreteAdaptTC::SetDiscreteTargetBase(const GridFunction &tspec_)
 {
    const int vdim = tspec_.FESpace()->GetVDim(),
              ndof = tspec_.FESpace()->GetNDofs();
-
    ncomp += vdim;
 
    // need to append data to tspec
@@ -1950,7 +1949,7 @@ void DiscreteAdaptTC::SetTspecAtIndex(int idx, const GridFunction &tspec_)
 {
    const int vdim = tspec_.FESpace()->GetVDim(),
              ndof = tspec_.FESpace()->GetNDofs();
-   MFEM_VERIFY(ndof == tspec.Size()/ncomp, "Inconsistency in SetTargetSpec.");
+   MFEM_VERIFY(ndof == tspec.Size()/ncomp, "Inconsistency in SetTspecAtIndex.");
 
    const auto tspec__d = tspec_.Read();
    auto tspec_d = tspec.ReadWrite();
@@ -2698,16 +2697,18 @@ void DiscreteAdaptTC::ComputeElementTargetsGradient(const IntegrationRule &ir,
    Jtrcomp.Clear();
 }
 
-void DiscreteAdaptTC:: UpdateGradientTargetSpecification(const Vector &x,
-                                                         const real_t dx,
-                                                         bool reuse_flag,
-                                                         int x_ordering)
+void DiscreteAdaptTC::
+UpdateGradientTargetSpecification(const Vector &x, real_t dx,
+                                  bool reuse_flag, int x_ordering)
 {
    if (reuse_flag && good_tspec_grad) { return; }
 
    const int dim = tspec_fesv->GetFE(0)->GetDim(),
              cnt = x.Size()/dim;
 
+   MFEM_VERIFY(tspec_fesv->GetVSize() / ncomp == cnt,
+               "FD with discrete adaptivity assume mesh_order = field_order.");
+
    tspec_pert1h.SetSize(x.Size()*ncomp);
 
    Vector TSpecTemp;
@@ -2734,16 +2735,18 @@ void DiscreteAdaptTC:: UpdateGradientTargetSpecification(const Vector &x,
 }
 
 void DiscreteAdaptTC::
-UpdateHessianTargetSpecification(const Vector &x,real_t dx,
+UpdateHessianTargetSpecification(const Vector &x, real_t dx,
                                  bool reuse_flag, int x_ordering)
 {
-
    if (reuse_flag && good_tspec_hess) { return; }
 
    const int dim    = tspec_fesv->GetFE(0)->GetDim(),
              cnt    = x.Size()/dim,
              totmix = 1+2*(dim-2);
 
+   MFEM_VERIFY(tspec_fesv->GetVSize() / ncomp == cnt,
+               "FD with discrete adaptivity assume mesh_order = field_order.");
+
    tspec_pert2h.SetSize(cnt*dim*ncomp);
    tspec_pertmix.SetSize(cnt*totmix*ncomp);
 
@@ -2909,6 +2912,11 @@ void TMOP_Integrator::EnableAdaptiveLimiting(const GridFunction &z0,
                                              Coefficient &coeff,
                                              AdaptivityEvaluator &ae)
 {
+   const char* gf_fe_name = z0.FESpace()->FEColl()->Name();
+   const char* mesh_fe_name =
+      z0.FESpace()->GetMesh()->GetNodalFESpace()->FEColl()->Name();
+   MFEM_VERIFY(strcmp(gf_fe_name, mesh_fe_name) == 0,
+               "Incompatible FE spaces for the adaptive limiting field.");
    adapt_lim_gf0 = &z0;
    delete adapt_lim_gf;
    adapt_lim_gf   = new GridFunction(z0);
@@ -2926,6 +2934,11 @@ void TMOP_Integrator::EnableAdaptiveLimiting(const ParGridFunction &z0,
                                              Coefficient &coeff,
                                              AdaptivityEvaluator &ae)
 {
+   const char* gf_fe_name = z0.FESpace()->FEColl()->Name();
+   const char* mesh_fe_name =
+      z0.FESpace()->GetMesh()->GetNodalFESpace()->FEColl()->Name();
+   MFEM_VERIFY(strcmp(gf_fe_name, mesh_fe_name) == 0,
+               "Incompatible FE spaces for the adaptive limiting field.");
    adapt_lim_gf0 = &z0;
    adapt_lim_pgf0 = &z0;
    delete adapt_lim_gf;

fem/tmop_tools.cpp

@@ -118,7 +118,7 @@ void AdvectorCG::ComputeAtNewPositionScalar(const Vector &new_nodes,
       h_min = std::min(h_min, m->GetElementSize(i));
    }
    real_t v_max = 0.0;
-   const int s = new_field.Size();
+   const int s  = u.Size()/m->Dimension();
 
    u.HostReadWrite();
    for (int i = 0; i < s; i++)
@@ -181,7 +181,7 @@ void AdvectorCG::ComputeAtNewPositionScalar(const Vector &new_nodes,
 
    // Trim the overshoots and undershoots.
    new_field.HostReadWrite();
-   for (int i = 0; i < s; i++)
+   for (int i = 0; i < new_field.Size(); i++)
    {
       if (new_field(i) < glob_minv) { new_field(i) = glob_minv; }
       if (new_field(i) > glob_maxv) { new_field(i) = glob_maxv; }
@@ -348,8 +348,10 @@ void InterpolatorFP::SetInitialField(const Vector &init_nodes,
 {
    nodes0 = init_nodes;
    Mesh *m = mesh;
+   FiniteElementSpace *f = fes;
 #ifdef MFEM_USE_MPI
    if (pmesh) { m = pmesh; }
+   if (pfes)  { f = pfes; }
 #endif
    m->SetNodes(nodes0);
 
@@ -363,14 +365,9 @@ void InterpolatorFP::SetInitialField(const Vector &init_nodes,
       delete finder;
    }
 
-   FiniteElementSpace *f = fes;
 #ifdef MFEM_USE_MPI
-   if (pfes)
-   {
-      f = pfes;
-      finder = new FindPointsGSLIB(pfes->GetComm());
-   }
-   else { finder = new FindPointsGSLIB(); }
+   if (pfes) { finder = new FindPointsGSLIB(pfes->GetComm()); }
+   else      { finder = new FindPointsGSLIB(); }
 #else
    finder = new FindPointsGSLIB();
 #endif
@@ -378,13 +375,74 @@ void InterpolatorFP::SetInitialField(const Vector &init_nodes,
 
    field0_gf.SetSpace(f);
    field0_gf = init_field;
+
+   // Check if the mesh nodes and the field nodes coincide.
+   const bool nodes_mismatch = init_nodes.Size() / m->Dimension() !=
+                               field0_gf.Size()  / f->GetVDim();
+   if (nodes_mismatch)
+   {
+      delete fes_field_nodes;
+      fes_field_nodes = new FiniteElementSpace(m, f->FEColl(), m->Dimension());
+   }
 }
 
 void InterpolatorFP::ComputeAtNewPosition(const Vector &new_nodes,
                                           Vector &new_field,
                                           int new_nodes_ordering)
 {
-   finder->Interpolate(new_nodes, field0_gf, new_field, new_nodes_ordering);
+   // Get physical node locations corresponding to field0_gf
+   if (fes_field_nodes)
+   {
+      Vector mapped_nodes;
+      GetFieldNodesPosition(new_nodes, mapped_nodes);
+      finder->Interpolate(mapped_nodes, field0_gf, new_field,
+                          fes_field_nodes->GetOrdering());
+   }
+   else
+   {
+      finder->Interpolate(new_nodes, field0_gf, new_field, new_nodes_ordering);
+   }
+}
+
+void InterpolatorFP::GetFieldNodesPosition(const Vector &mesh_nodes,
+                                           Vector &nodes_pos) const
+{
+   MFEM_VERIFY(fes_field_nodes, "InterpolatorFP: fes_field_nodes is not set.");
+
+   Mesh *m = fes_field_nodes->GetMesh();
+   const int nelem     = fes_field_nodes->GetNE();
+   const int n_f_nodes = fes_field_nodes->GetNDofs();
+   const int dim       = m->Dimension();
+   if (nelem == 0) { return; }
+   Array<int> dofs;
+   Vector e_xyz;
+   nodes_pos.SetSize(n_f_nodes*dim);
+   const FiniteElementSpace *mesh_fes = m->GetNodalFESpace();
+
+   for (int e = 0; e < nelem; e++)
+   {
+      mesh_fes->GetElementVDofs(e, dofs);
+      int n_mdofs = dofs.Size()/dim;
+      mesh_nodes.GetSubVector(dofs, e_xyz); //e_xyz is ordered by nodes here
+      const FiniteElement *mfe = mesh_fes->GetFE(e);
+      Vector shape(n_mdofs);
+
+      auto ir = fes_field_nodes->GetFE(e)->GetNodes();
+      const int n_gf_pts = ir.GetNPoints();
+      Vector gf_xyz(n_gf_pts*dim);
+      for (int q = 0; q < n_gf_pts; q++)
+      {
+         IntegrationPoint ip = ir.IntPoint(q);
+         mfe->CalcShape(ip, shape);
+         for (int d = 0; d < dim; d++)
+         {
+            Vector x(e_xyz.GetData() + d*n_mdofs, n_mdofs);
+            gf_xyz(d*n_gf_pts + q) = x*shape; // order by nodes
+         }
+      }
+      fes_field_nodes->GetElementVDofs(e, dofs);
+      nodes_pos.SetSubVector(dofs, gf_xyz);
+   }
 }
 
 #endif
@@ -437,7 +495,7 @@ real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
       // Check for convergence
       if (init_fit_max_err < surf_fit_max_err_limit)
       {
-         if (print_options.iterations)
+         if (print_options.iterations || print_options.warnings)
          {
             mfem::out << "TMOPNewtonSolver converged "
                       "based on the surface fitting error.\n";

fem/tmop_tools.hpp

@@ -58,8 +58,14 @@ private:
    Vector nodes0;
    GridFunction field0_gf;
    FindPointsGSLIB *finder;
+   // FE space for the nodes of the solution GridFunction.
+   FiniteElementSpace *fes_field_nodes;
+
+   void GetFieldNodesPosition(const Vector &mesh_nodes,
+                              Vector &nodes_pos) const;
+
 public:
-   InterpolatorFP() : finder(NULL) { }
+   InterpolatorFP() : finder(NULL), fes_field_nodes(NULL) { }
 
    void SetInitialField(const Vector &init_nodes,
                         const Vector &init_field) override;
@@ -77,6 +83,7 @@ public:
    {
       finder->FreeData();
       delete finder;
+      delete fes_field_nodes;
    }
 };
 #endif

fem/transfer.cpp

@@ -1159,10 +1159,10 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space()
    }
 
    // need scalar to keep dimensions matching (operators are built to apply individually on each vdim)
-   fes_ho_scalar = new FiniteElementSpace(fes_ho.GetMesh(),
-                                          fes_ho.FEColl(), 1);
-   fes_lor_scalar = new FiniteElementSpace(fes_lor.GetMesh(),
-                                           fes_lor.FEColl(), 1);
+   fes_ho_scalar.reset(new FiniteElementSpace(fes_ho.GetMesh(),
+                                              fes_ho.FEColl(), 1));
+   fes_lor_scalar.reset(new FiniteElementSpace(fes_lor.GetMesh(),
+                                               fes_lor.FEColl(), 1));
 
    const CoarseFineTransformations& cf_tr = mesh_lor->GetRefinementTransforms();
 
@@ -1197,7 +1197,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space()
    ML_inv_ea.SetSize(ndof_lor);
    ML_inv_ea = 0.0;
 
-   BilinearForm Mho(fes_ho_scalar);
+   BilinearForm Mho(fes_ho_scalar.get());
    Mho.SetAssemblyLevel(AssemblyLevel::PARTIAL);
    Mho.AddDomainIntegrator(new MassIntegrator);
    Mho.Assemble();
@@ -1207,7 +1207,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space()
    M_H = 0.0;
    Mho.Mult(ones_ho, M_H);
 
-   BilinearForm Mlor(fes_lor_scalar);
+   BilinearForm Mlor(fes_lor_scalar.get());
    Mlor.SetAssemblyLevel(AssemblyLevel::PARTIAL);
    Mlor.AddDomainIntegrator(new MassIntegrator);
    Mlor.Assemble();
@@ -1236,7 +1236,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space()
    precon.reset(new OperatorJacobiSmoother(M_H, ess_tdof_list));
 
    TransposeOperator* RT = new TransposeOperator(R.get());
-   RTxM_LH.reset(new ProductOperator(RT, M_LH.get(), false, false));
+   RTxM_LH.reset(new ProductOperator(RT, M_LH.get(), true, false));
 
    SetupPCG();
 }
@@ -1260,10 +1260,10 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space
    }
 
    // need scalar to keep dimensions matching (operators are built to apply individually on each vdim)
-   pfes_ho_scalar = new ParFiniteElementSpace(pfes_ho.GetParMesh(),
-                                              pfes_ho.FEColl(), 1);
-   pfes_lor_scalar = new ParFiniteElementSpace(pfes_lor.GetParMesh(),
-                                               pfes_lor.FEColl(), 1);
+   pfes_ho_scalar.reset(new ParFiniteElementSpace(pfes_ho.GetParMesh(),
+                                                  pfes_ho.FEColl(), 1));
+   pfes_lor_scalar.reset(new ParFiniteElementSpace(pfes_lor.GetParMesh(),
+                                                   pfes_lor.FEColl(), 1));
 
    const CoarseFineTransformations& cf_tr = mesh_lor->GetRefinementTransforms();
 
@@ -1292,7 +1292,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space
    ML_inv_ea.SetSize(ndof_lor);
    ML_inv_ea = 0.0;
 
-   BilinearForm pMho(pfes_ho_scalar);
+   ParBilinearForm pMho(pfes_ho_scalar.get());
    pMho.SetAssemblyLevel(AssemblyLevel::PARTIAL);
    pMho.AddDomainIntegrator(new MassIntegrator);
    pMho.Assemble();
@@ -1302,7 +1302,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space
    M_H = 0.0;
    pMho.Mult(ones_ho, M_H);
 
-   ParBilinearForm pMlor(pfes_lor_scalar);
+   ParBilinearForm pMlor(pfes_lor_scalar.get());
    pMlor.SetAssemblyLevel(AssemblyLevel::PARTIAL);
    pMlor.AddDomainIntegrator(new MassIntegrator);
    pMlor.Assemble();
@@ -1318,10 +1318,11 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space
    // **************************
    // mixed mass M_LH
    // **************************
-   MixedMassEA(*pfes_ho_scalar, *pfes_lor_scalar, M_LH_ea, d_mt);
+   MixedMassEA(*pfes_ho_scalar.get(), *pfes_lor_scalar.get(), M_LH_ea, d_mt);
 
    // Set ownership
-   M_LH_local_op = new H1SpaceMixedMassOperator(pfes_ho_scalar, pfes_lor_scalar,
+   M_LH_local_op = new H1SpaceMixedMassOperator(pfes_ho_scalar.get(),
+                                                pfes_lor_scalar.get(),
                                                 &ho2lor, &M_LH_ea);
 
    const Operator *P_ho = pfes_ho_scalar->GetProlongationMatrix();
@@ -1336,10 +1337,11 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space
          Operator *Pt_lor = new TransposeOperator(P_lor);
          RML_inv.SetSize(pfes_lor_scalar->GetTrueVSize());
          GetTDofs(*pfes_lor_scalar, ML_inv_ea, RML_inv);
-         ML_inv_vea.reset(new H1SpaceLumpedMassOperator(pfes_ho_scalar, pfes_lor_scalar,
+         ML_inv_vea.reset(new H1SpaceLumpedMassOperator(pfes_ho_scalar.get(),
+                                                        pfes_lor_scalar.get(),
                                                         RML_inv));
-         M_LH.reset(new TripleProductOperator(Pt_lor, M_LH_local_op, P_ho, false,
-                                              false, false));
+         M_LH.reset(new TripleProductOperator(Pt_lor, M_LH_local_op, P_ho, true,
+                                              true, false));
 
          Vector RM_H(pfes_ho_scalar->GetTrueVSize());
          GetTDofsTranspose(*pfes_ho_scalar, M_H, RM_H);
@@ -1347,12 +1349,13 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space
       }
       else if (P_ho)
       {
-         ML_inv_vea.reset(new H1SpaceLumpedMassOperator(pfes_ho_scalar, pfes_lor_scalar,
+         ML_inv_vea.reset(new H1SpaceLumpedMassOperator(pfes_ho_scalar.get(),
+                                                        pfes_lor_scalar.get(),
                                                         ML_inv_ea));
-         M_LH.reset(new ProductOperator(M_LH_local_op, P_ho, false, false));
+         M_LH.reset(new ProductOperator(M_LH_local_op, P_ho, true, false));
 
          Vector RM_H(pfes_ho_scalar->GetTrueVSize());
-         GetTDofsTranspose(*pfes_ho_scalar, M_H, RM_H);
+         GetTDofsTranspose(*pfes_ho_scalar.get(), M_H, RM_H);
          precon.reset(new OperatorJacobiSmoother(RM_H, ess_tdof_list));
       }
       else if (P_lor)
@@ -1360,9 +1363,10 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space
          Operator *Pt_lor = new TransposeOperator(P_lor);
          RML_inv.SetSize(pfes_lor_scalar->GetTrueVSize());
          GetTDofsTranspose(*pfes_lor_scalar, ML_inv_ea, RML_inv);
-         ML_inv_vea.reset(new H1SpaceLumpedMassOperator(pfes_ho_scalar, pfes_lor_scalar,
+         ML_inv_vea.reset(new H1SpaceLumpedMassOperator(pfes_ho_scalar.get(),
+                                                        pfes_lor_scalar.get(),
                                                         RML_inv));
-         M_LH.reset(new ProductOperator(Pt_lor, M_LH_local_op, false, false));
+         M_LH.reset(new ProductOperator(Pt_lor, M_LH_local_op, true, true));
          R.reset(new ProductOperator(ML_inv_vea.get(), M_LH.get(), false,
                                      false));
 
@@ -1370,7 +1374,8 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space
       }
       else
       {
-         ML_inv_vea.reset(new H1SpaceLumpedMassOperator(pfes_ho_scalar, pfes_lor_scalar,
+         ML_inv_vea.reset(new H1SpaceLumpedMassOperator(pfes_ho_scalar.get(),
+                                                        pfes_lor_scalar.get(),
                                                         ML_inv_ea));
          M_LH.reset(M_LH_local_op);
 
@@ -1381,7 +1386,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space
                                false));
 
    TransposeOperator* RT = new TransposeOperator(R.get());
-   RTxM_LH.reset(new ProductOperator(RT, M_LH.get(), false, false));
+   RTxM_LH.reset(new ProductOperator(RT, M_LH.get(), true, false));
 
    SetupPCG();
 }

fem/transfer.hpp

@@ -472,16 +472,16 @@ public:
       Operator *M_LH_local_op;
 
       // Scalar finite element spaces for stored Tdof-to-and-from-LDof maps.
-      FiniteElementSpace* fes_ho_scalar;
-      FiniteElementSpace* fes_lor_scalar;
+      std::unique_ptr<FiniteElementSpace> fes_ho_scalar;
+      std::unique_ptr<FiniteElementSpace> fes_lor_scalar;
       // Element Assembled mixed mass
       Vector M_LH_ea;
       // Element Assembled lumped M_L inverse built via EA. Stores diagonal as a Ldof vector.
       Vector ML_inv_ea;
 
 #ifdef MFEM_USE_MPI
-      ParFiniteElementSpace* pfes_ho_scalar;
-      ParFiniteElementSpace* pfes_lor_scalar;
+      std::unique_ptr<ParFiniteElementSpace> pfes_ho_scalar;
+      std::unique_ptr<ParFiniteElementSpace> pfes_lor_scalar;
       Vector RML_inv;
 #endif
 

general/CMakeLists.txt

@@ -45,6 +45,7 @@ list(APPEND HDRS
   gecko.hpp
   globals.hpp
   zstr.hpp
+  handle.hpp
   hash.hpp
   isockstream.hpp
   kdtree.hpp

general/array.cpp

@@ -112,7 +112,7 @@ void Array<T>::PartialSum()
 
 // Sum
 template <class T>
-T Array<T>::Sum()
+T Array<T>::Sum() const
 {
    T sum = static_cast<T>(0);
    for (int i = 0; i < size; i++)

general/array.hpp

@@ -23,6 +23,7 @@
 #include <cstring>
 #include <algorithm>
 #include <type_traits>
+#include <initializer_list>
 
 namespace mfem
 {
@@ -88,24 +89,14 @@ public:
    template <typename CT>
    inline Array(const Array<CT> &src);
 
-   /// Deep copy from a braced init-list of convertible type
+   /// Construct an Array from a C-style array of static length
    template <typename CT, int N>
    explicit inline Array(const CT (&values)[N]);
 
-   /**
-    * @brief Construct a new Array object from an initializer list.
-    *
-    * @param init_list List of entities to construct from.
-    */
-   Array(const std::initializer_list<T> &init_list)
-      : Array(static_cast<int>(init_list.size()))
-   {
-      auto * it = GetData();
-      for (auto value : init_list)
-      {
-         *it++ = value;
-      }
-   }
+   /// Construct an Array from a braced initializer list of convertible type
+   template <typename CT, typename std::enable_if<
+                std::is_convertible<CT,T>::value,bool>::type = true>
+   explicit inline Array(std::initializer_list<CT> values);
 
    /// Move constructor ("steals" data from 'src')
    inline Array(Array<T> &&src) { Swap(src, *this); }
@@ -295,7 +286,7 @@ public:
    void Unique()
    {
       T* end = std::unique((T*)data, data + size);
-      SetSize(end - data);
+      SetSize((int)(end - data));
    }
 
    /// Return 1 if the array is sorted from lowest to highest.  Otherwise return 0.
@@ -308,7 +299,7 @@ public:
    void PartialSum();
 
    /// Return the sum of all the array entries using the '+'' operator for class 'T'.
-   T Sum();
+   T Sum() const;
 
    /// Set all entries of the array to the provided constant.
    inline void operator=(const T &a);
@@ -701,10 +692,18 @@ inline Array<T>::Array(const Array<CT> &src)
    for (int i = 0; i < size; i++) { (*this)[i] = T(src[i]); }
 }
 
+template <typename T>
+template <typename CT, typename std::enable_if<
+             std::is_convertible<CT,T>::value,bool>::type>
+inline Array<T>::Array(std::initializer_list<CT> values) : Array(values.size())
+{
+   std::copy(values.begin(), values.end(), begin());
+}
+
 template <typename T> template <typename CT, int N>
 inline Array<T>::Array(const CT (&values)[N]) : Array(N)
 {
-   for (int i = 0; i < size; i++) { (*this)[i] = T(values[i]); }
+   std::copy(values, values + N, begin());
 }
 
 template <class T>
@@ -902,7 +901,7 @@ inline int Array<T>::FindSorted(const T &el) const
    const T *begin = data, *end = begin + size;
    const T* first = std::lower_bound(begin, end, el);
    if (first == end || !(*first == el)) { return  -1; }
-   return first - begin;
+   return (int)(first - begin);
 }
 
 template <class T>

general/arrays_by_name.hpp

@@ -64,7 +64,7 @@ public:
    ArraysByName(ArraysByName &&src) noexcept = default;
 
    /// Return the number of named arrays in the container
-   int Size() const { return data.size(); }
+   int Size() const { return static_cast<int>(data.size()); }
 
    /// Return an STL set of strings giving the names of the arrays
    inline std::set<std::string> GetNames() const;

general/communication.hpp

@@ -431,8 +431,8 @@ struct VarMessage
    void Isend(int rank, MPI_Comm comm)
    {
       Encode(rank);
-      MPI_Isend((void*) data.data(), data.length(), MPI_BYTE, rank, Tag, comm,
-                &send_request);
+      MPI_Isend((void*) data.data(), static_cast<int>(data.length()), MPI_BYTE, rank,
+                Tag, comm, &send_request);
    }
 
    /** @brief Non-blocking synchronous send to processor 'rank'.
@@ -441,8 +441,8 @@ struct VarMessage
    void Issend(int rank, MPI_Comm comm)
    {
       Encode(rank);
-      MPI_Issend((void*) data.data(), data.length(), MPI_BYTE, rank, Tag, comm,
-                 &send_request);
+      MPI_Issend((void*) data.data(), static_cast<int>(data.length()), MPI_BYTE, rank,
+                 Tag, comm, &send_request);
    }
 
    /// Helper to send all messages in a rank-to-message map container.
@@ -538,7 +538,7 @@ struct VarMessage
    template<typename MapT>
    static void RecvAll(MapT& rank_msg, MPI_Comm comm)
    {
-      int recv_left = rank_msg.size();
+      int recv_left = static_cast<int>(rank_msg.size());
       while (recv_left > 0)
       {
          int rank, size;

general/device.cpp

@@ -150,6 +150,9 @@ Device::Device()
 
 Device::~Device()
 {
+#ifdef MFEM_USE_MPI
+   Hypre::Finalize();
+#endif
    if ( device_env && !destroy_mm) { return; }
    if (!device_env &&  destroy_mm && !mem_host_env)
    {
@@ -255,7 +258,15 @@ void Device::Configure(const std::string &device, const int device_id)
    destroy_mm = true;
 
 #ifdef MFEM_USE_MPI
-   Hypre::InitDevice();
+#if defined(HYPRE_USING_GPU) && (MFEM_HYPRE_VERSION >= 23100)
+   // Skip the call to Hypre::InitDevice() if HYPRE is not initialized, e.g.
+   // * if running a serial code
+   // * if running with the environment variable MFEM_DEVICE set.
+   if (HYPRE_Initialized())
+   {
+      Hypre::InitDevice();
+   }
+#endif
 #endif
 }
 

general/globals.cpp

@@ -9,12 +9,17 @@
 // terms of the BSD-3 license. We welcome feedback and contributions, see file
 // CONTRIBUTING.md for details.
 
+#ifdef _WIN32
+// Turn off CRT deprecation warnings for getenv
+#define _CRT_SECURE_NO_WARNINGS
+#endif
 
 #include "../config/config.hpp"
 #include "globals.hpp"
 #include <iostream>
 #include <sstream>
 #include <iomanip>
+#include <cstdlib>  // getenv
 
 namespace mfem
 {
@@ -71,4 +76,8 @@ void SetGlobalMPI_Comm(MPI_Comm comm)
 
 #endif
 
+const char* getenv(const char* name)
+{
+   return ::getenv(name);
+}
 }

general/globals.hpp

@@ -107,6 +107,8 @@ void SetGlobalMPI_Comm(MPI_Comm comm);
 
 #endif
 
+const char* getenv(const char* name);
+
 } // namespace mfem
 
 #endif

general/handle.hpp

@@ -0,0 +1,200 @@
+// Copyright (c) 2010-2024, 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_HANDLE_HPP
+#define MFEM_HANDLE_HPP
+
+#include "../config/config.hpp"
+#include <memory>
+
+namespace mfem
+{
+
+/// @brief A smart pointer class that may represent either shared ownership, or
+/// a non-owning borrow.
+///
+/// A Handle may either be owning or non-owning. Non-owning Handle%s point to
+/// externally owned data; it is the responsibility of the user both to ensure
+/// that the data remains valid as long as the Handle is alive and to delete the
+/// pointer when it is no longer needed. Owning Handle%s use <a
+/// href="https://en.cppreference.com/w/cpp/memory/shared_ptr">
+/// std::shared_ptr</a> to implement reference counting. The underlying data
+/// will be valid as long as there is at least one live copy. When the last
+/// Handle is destroyed, the pointer is deleted.
+///
+/// Both types of Handle%s can be copied, moved, stored in standard containers,
+/// etc.
+///
+/// A non-owning Handle may assume ownership over its data, but an owning Handle
+/// cannot release ownership over its data.
+///
+/// It is an invariant of this class that **at most** one of the data members
+/// @a not_owned and @a owned will be non-null.
+template <typename T>
+class Handle
+{
+   /// If this is a non-owning handle, @a not_owned will point to the data.
+   T *not_owned = nullptr;
+
+   /// If this is an owning handle, @a owned will point to the data.
+   std::shared_ptr<T> owned = nullptr;
+
+   /// @brief Types @a Handle<T> and @a %Handle\<U\> are friends to allow
+   /// construction of one from another when @a T and @a U are convertible
+   /// types.
+   template <typename U> friend class Handle;
+
+public:
+   /// Create an empty (null) Handle.
+   Handle() = default;
+
+   /// @brief Create a Handle pointing to @a t.
+   ///
+   /// If @a take_ownership is true, then the Handle assumes ownership over the
+   /// pointer, and it should not be deleted externally. Otherwise, the Handle
+   /// will be non-owning, and it is the user's responsibility to ensure the
+   /// correct lifetime of @a t.
+   Handle(T *t, bool take_ownership)
+   {
+      if (take_ownership) { owned.reset(t); }
+      else { not_owned = t; }
+   }
+
+   /// Create a Handle from a std::shared_ptr (sharing ownership with @a t).
+   Handle(const std::shared_ptr<T> &t) : owned(t) { }
+
+   /// @brief Copy constructor.
+   ///
+   /// Copying an owning Handle results in another owning handle. Copying a
+   /// non-owning handle results in a non-owning handle.
+   Handle(const Handle &other) = default;
+
+   /// Move constructor (see Handle(const Handle&)).
+   Handle(Handle &&other) = default;
+
+   /// @brief Constructs a copy of @a u, where type @a U is convertible to @a T.
+   ///
+   /// This allows the construction of Handle<Base> from Handle<Derived>.
+   template <typename U>
+   Handle(const Handle<U> &u) : not_owned(u.not_owned), owned(u.owned) { }
+
+   /// @brief Move-constructs from @a u, where type @a U is convertible to @a T.
+   ///
+   /// See @ref Handle(const Handle<U>&).
+   template <typename U>
+   Handle(Handle<U> &&u) : not_owned(u.not_owned), owned(u.owned) { }
+
+   /// Destructor. If the Handle is owning, decrement the reference count.
+   ~Handle() = default;
+
+   /// Copy assignment (see Handle(const Handle&)).
+   Handle &operator=(const Handle &other) = default;
+
+   /// Move assignment (see Handle(const Handle&)).
+   Handle &operator=(Handle &&other) = default;
+
+   /// Returns the contained pointer (may be null).
+   T *Get() const
+   {
+      if (not_owned) { return not_owned; }
+      else { return owned.get(); }
+   }
+
+   /// @brief If the Handle is owning, return a copy of the underlying shared
+   /// pointer.
+   ///
+   /// @warning If the Handle is non-owning (even if non-null), this will return
+   /// and empty (null) shared pointer.
+   std::shared_ptr<T> GetSharedPtr() const { return owned; }
+
+   /// Dereference operator. The Handle must be non-null.
+   T &operator*() const { return *Get(); }
+
+   /// Member access (arrow) operator. The Handle must be non-null.
+   T *operator->() const { return Get(); }
+
+   /// @brief Returns true if the Handle is owning, false if it is non-owning.
+   ///
+   /// Returns false if the Handle is null (empty).
+   bool IsOwner() const { return owned; }
+
+   /// Returns true if the Handle is non-null.
+   explicit operator bool() const { return not_owned || owned; }
+
+   /// @brief Assume owernship of the data.
+   ///
+   /// If the Handle is already owning, this does nothing.
+   void MakeOwner()
+   {
+      if (owned) { return; }
+      owned.reset(not_owned);
+      not_owned = nullptr;
+   }
+
+   /// @brief Reset the Handle to be empty.
+   ///
+   /// If the Handle is owning, this will decrement the reference count.
+   void Reset()
+   {
+      owned.reset();
+      not_owned = nullptr;
+   }
+
+   /// @brief Reset the Handle to point to @a t.
+   ///
+   /// The Handle may assume ownership of the pointer according to @a
+   /// take_ownership (see @ref Handle(T*, bool)).
+   void Reset(T *t, bool take_ownership)
+   {
+      if (take_ownership)
+      {
+         owned.reset(t);
+         not_owned = nullptr;
+      }
+      else
+      {
+         owned.reset();
+         not_owned = t;
+      }
+   }
+
+   /// Reset the Handle to share ownership with @a t.
+   void Reset(const std::shared_ptr<T> &t)
+   {
+      owned = t;
+      not_owned = nullptr;
+   }
+};
+
+/// @brief Return a new owning Handle, where the pointed-to object is a new
+/// object constructed using the given arguments.
+///
+/// This is analogous to <a
+/// href="https://en.cppreference.com/w/cpp/memory/shared_ptr/make_shared">
+/// std::make_shared</a>.
+template <typename T, typename... Args>
+Handle<T> MakeOwning(Args&&... args)
+{
+   T *t = new T(std::forward<Args>(args)...);
+   return Handle<T>(t, true);
+}
+
+/// Return a new owning Handle pointing to @a t.
+template <typename T>
+Handle<T> Owning(T *t) { return Handle<T>(t, true); }
+
+/// Return a new non-owning Handle pointing to @a t.
+template <typename T>
+Handle<T> NonOwning(T *t) { return Handle<T>(t, false); }
+
+} // namespace mfem
+
+#endif

general/isockstream.cpp

@@ -87,7 +87,7 @@ int isockstream::establish()
 #if defined(__APPLE__)
       if (bind(sfd, (const struct sockaddr *)rp->ai_addr, rp->ai_addrlen) < 0)
 #else
-      if (bind(sfd, rp->ai_addr, rp->ai_addrlen) < 0)
+      if (bind(sfd, rp->ai_addr, static_cast<socklen_t>(rp->ai_addrlen)) < 0)
 #endif
       {
          mfem::err << "isockstream::establish(): bind() failed!" << std::endl;

general/socketstream.cpp

@@ -10,7 +10,7 @@
 // CONTRIBUTING.md for details.
 
 #ifdef _WIN32
-// Turn off CRT deprecation warnings for strerror (VS 2013)
+// Turn off CRT deprecation warnings for strerror
 #define _CRT_SECURE_NO_WARNINGS
 #endif
 
@@ -30,13 +30,13 @@
 #include <ws2tcpip.h>
 #ifdef _MSC_VER
 typedef int ssize_t;
+typedef int socklen_t;
 // Link with ws2_32.lib
 #pragma comment(lib, "ws2_32.lib")
 #endif
 #endif
 
 #ifdef MFEM_USE_GNUTLS
-#include <cstdlib>  // getenv
 #ifndef MFEM_USE_GNUTLS_X509
 #include <gnutls/openpgp.h>
 #endif
@@ -146,7 +146,8 @@ int socketbuf::open(const char hostname[], int port)
       }
 #endif
 
-      if (connect(socket_descriptor, rp->ai_addr, rp->ai_addrlen) < 0)
+      if (connect(socket_descriptor, rp->ai_addr,
+                  static_cast<socklen_t>(rp->ai_addrlen)) < 0)
       {
          closesocket(socket_descriptor);
          socket_descriptor = -2;

linalg/CMakeLists.txt

@@ -24,7 +24,7 @@ list(APPEND SRCS
   constraints.cpp
   densemat.cpp
   symmat.cpp
-  handle.cpp
+  op_handle.cpp
   matrix.cpp
   ode.cpp
   operator.cpp
@@ -51,7 +51,7 @@ list(APPEND HDRS
   dinvariants.hpp
   symmat.hpp
   dtensor.hpp
-  handle.hpp
+  op_handle.hpp
   invariants.hpp
   kernels.hpp
   lapack.hpp

linalg/hypre.cpp

@@ -28,16 +28,23 @@ namespace mfem
 {
 
 bool Hypre::configure_runtime_policy_from_mfem = true;
+Hypre::State Hypre::state = Hypre::State::UNINITIALIZED;
 
-Hypre::Hypre()
+void Hypre::Init()
 {
+   if (state != State::INITIALIZED)
+   {
 #if MFEM_HYPRE_VERSION >= 21900
-   // Initializing hypre
-   HYPRE_Init();
+      HYPRE_Init();
 #endif
-
-   // Global hypre options that we set by default
-   SetDefaultOptions();
+      SetDefaultOptions();
+      // Apply the setting of 'configure_runtime_policy_from_mfem' according to
+      // the current configuration of the mfem::Device (HYPRE >= 2.31.0):
+      InitDevice();
+      // Create the singleton Hypre object AFTER initializing HYPRE:
+      Instance();
+   }
+   state = State::INITIALIZED;
 }
 
 void Hypre::InitDevice()
@@ -48,6 +55,8 @@ void Hypre::InitDevice()
 #if defined(HYPRE_USING_GPU) && (MFEM_HYPRE_VERSION >= 23100)
    if (configure_runtime_policy_from_mfem)
    {
+      MFEM_VERIFY(HYPRE_Initialized(), "HYPRE must be initialized before"
+                  " calling Hypre::InitDevice()");
       if (Device::Allows(Backend::DEVICE_MASK & ~Backend::DEBUG_DEVICE))
       {
          HYPRE_SetMemoryLocation(HYPRE_MEMORY_DEVICE);
@@ -65,14 +74,13 @@ void Hypre::InitDevice()
 
 void Hypre::Finalize()
 {
-   Hypre &hypre = Instance();
-   if (!hypre.finalized)
+   if (state != State::UNINITIALIZED)
    {
 #if MFEM_HYPRE_VERSION >= 21900
       HYPRE_Finalize();
 #endif
-      hypre.finalized = true;
    }
+   state = State::UNINITIALIZED;
 }
 
 void Hypre::SetDefaultOptions()
@@ -5355,8 +5363,8 @@ void HypreBoomerAMG::SetAdvectiveOptions(int distanceR,
    int ns_down = 0, ns_up = 0, ns_coarse; // init to suppress gcc warnings
    if (distanceR > 0)
    {
-      ns_down = prerelax.length();
-      ns_up = postrelax.length();
+      ns_down = static_cast<int>(prerelax.length());
+      ns_up = static_cast<int>(postrelax.length());
       ns_coarse = 1;
 
       // Array to store relaxation scheme and pass to Hypre

linalg/hypre.hpp

@@ -68,10 +68,11 @@ class Hypre
 public:
    /// @brief Initialize hypre by calling HYPRE_Init() and set default options.
    /// After calling Hypre::Init(), hypre will be finalized automatically at
-   /// program exit.
+   /// program exit. May be re-initialized after finalize.
    ///
-   /// Calling HYPRE_Finalize() manually is not compatible with this class.
-   static void Init() { Instance(); }
+   /// Calling HYPRE_Init() or HYPRE_Finalize() manually is only supported for
+   /// HYPRE 2.29.0+
+   static void Init();
 
    /// @brief Configure HYPRE's compute and memory policy.
    ///
@@ -94,6 +95,9 @@ public:
    ///
    /// Multiple calls to Hypre::Finalize() have no effect. This function can be
    /// called manually to more precisely control when hypre is finalized.
+   ///
+   /// Calling HYPRE_Init() or HYPRE_Finalize() manually is only supported for
+   /// HYPRE 2.29.0+
    static void Finalize();
 
    /// @brief Use MFEM's device policy to configure HYPRE's device policy, true
@@ -104,14 +108,20 @@ public:
    static bool configure_runtime_policy_from_mfem;
 
 private:
-   /// Calls HYPRE_Init() when the singleton is constructed.
-   Hypre();
+   /// Default constructor. Singleton object; private.
+   Hypre() = default;
+
+   /// Copy constructor. Deleted.
+   Hypre(Hypre&) = delete;
+
+   /// Move constructor. Deleted.
+   Hypre(Hypre&&) = delete;
 
    /// The singleton destructor (called at program exit) finalizes hypre.
    ~Hypre() { Finalize(); }
 
    /// Set the default hypre global options (mostly GPU-relevant).
-   void SetDefaultOptions();
+   static void SetDefaultOptions();
 
    /// Create and return the Hypre singleton object.
    static Hypre &Instance()
@@ -120,7 +130,10 @@ private:
       return hypre;
    }
 
-   bool finalized = false; ///< Has Hypre::Finalize() been called already?
+   enum class State { UNINITIALIZED, INITIALIZED };
+
+   /// Tracks whether Hypre was initialized or finalized by this class.
+   static State state;
 };
 
 

linalg/linalg.hpp

@@ -28,7 +28,7 @@
 #include "symmat.hpp"
 #include "ode.hpp"
 #include "solvers.hpp"
-#include "handle.hpp"
+#include "op_handle.hpp"
 #include "invariants.hpp"
 #include "constraints.hpp"
 #include "auxiliary.hpp"

linalg/op_handle.cpp

@@ -9,7 +9,7 @@
 // terms of the BSD-3 license. We welcome feedback and contributions, see file
 // CONTRIBUTING.md for details.
 
-#include "handle.hpp"
+#include "op_handle.hpp"
 #include "sparsemat.hpp"
 #ifdef MFEM_USE_MPI
 #include "petsc.hpp"

linalg/op_handle.hpp

@@ -9,8 +9,8 @@
 // terms of the BSD-3 license. We welcome feedback and contributions, see file
 // CONTRIBUTING.md for details.
 
-#ifndef MFEM_HANDLE_HPP
-#define MFEM_HANDLE_HPP
+#ifndef MFEM_OP_HANDLE_HPP
+#define MFEM_OP_HANDLE_HPP
 
 #include "../config/config.hpp"
 #include "operator.hpp"

linalg/operator.cpp

@@ -183,8 +183,8 @@ Operator * Operator::SetupRAP(const Operator *Pi, const Operator *Po)
    {
       if (!IsIdentityProlongation(Po))
       {
-         TransposeOperator * PoT = new TransposeOperator(Po);
-         rap = new ProductOperator(PoT, this, true,false);
+         rap = new ProductOperator(Owning(new TransposeOperator(Po)),
+                                   NonOwning(this));
       }
       else
       {
@@ -365,11 +365,10 @@ void SecondOrderTimeDependentOperator::ImplicitSolve(const real_t dt0,
    mfem_error("SecondOrderTimeDependentOperator::ImplicitSolve() is not overridden!");
 }
 
-SumOperator::SumOperator(const Operator *A, const real_t alpha,
-                         const Operator *B, const real_t beta,
-                         bool ownA, bool ownB)
-   : Operator(A->Height(), A->Width()),
-     A(A), B(B), alpha(alpha), beta(beta), ownA(ownA), ownB(ownB),
+SumOperator::SumOperator(Handle<const Operator> A_, const real_t alpha,
+                         Handle<const Operator> B_, const real_t beta)
+   : Operator(A_->Height(), A_->Width()),
+     A(A_), B(B_), alpha(alpha), beta(beta),
      z(A->Height())
 {
    MFEM_VERIFY(A->Width() == B->Width(),
@@ -381,53 +380,43 @@ SumOperator::SumOperator(const Operator *A, const real_t alpha,
                << "A->Height() = " << A->Height()
                << ", B->Height() = " << B->Height() );
 
+   if (auto SolverA = dynamic_cast<const Solver*>(A.Get()))
    {
-      const Solver* SolverA = dynamic_cast<const Solver*>(A);
-      const Solver* SolverB = dynamic_cast<const Solver*>(B);
-      if (SolverA)
-      {
-         MFEM_VERIFY(!(SolverA->iterative_mode),
-                     "Operator A of a SumOperator should not be in iterative mode");
-      }
-      if (SolverB)
-      {
-         MFEM_VERIFY(!(SolverB->iterative_mode),
-                     "Operator B of a SumOperator should not be in iterative mode");
-      }
+      MFEM_VERIFY(!(SolverA->iterative_mode),
+                  "Operator A of a SumOperator should not be in iterative mode");
+   }
+   if (auto SolverB = dynamic_cast<const Solver*>(B.Get()))
+   {
+      MFEM_VERIFY(!(SolverB->iterative_mode),
+                  "Operator B of a SumOperator should not be in iterative mode");
    }
-
 }
 
-SumOperator::~SumOperator()
-{
-   if (ownA) { delete A; }
-   if (ownB) { delete B; }
-}
+SumOperator::SumOperator(const Operator *A_, const real_t alpha,
+                         const Operator *B_, const real_t beta,
+                         bool own_A, bool own_B)
+   : SumOperator({A_, own_A}, alpha, {B_, own_B}, beta) { }
 
-ProductOperator::ProductOperator(const Operator *A, const Operator *B,
-                                 bool ownA, bool ownB)
-   : Operator(A->Height(), B->Width()),
-     A(A), B(B), ownA(ownA), ownB(ownB), z(A->Width())
+
+ProductOperator::ProductOperator(Handle<const Operator> A_,
+                                 Handle<const Operator> B_)
+   : Operator(A_->Height(), B_->Width()),
+     A(A_), B(B_), z(A->Width())
 {
    MFEM_VERIFY(A->Width() == B->Height(),
                "incompatible Operators: A->Width() = " << A->Width()
                << ", B->Height() = " << B->Height());
 
+   if (auto SolverB = dynamic_cast<const Solver*>(B.Get()))
    {
-      const Solver* SolverB = dynamic_cast<const Solver*>(B);
-      if (SolverB)
-      {
-         MFEM_VERIFY(!(SolverB->iterative_mode),
-                     "Operator B of a ProductOperator should not be in iterative mode");
-      }
+      MFEM_VERIFY(!(SolverB->iterative_mode),
+                  "Operator B of a ProductOperator should not be in iterative mode");
    }
 }
 
-ProductOperator::~ProductOperator()
-{
-   if (ownA) { delete A; }
-   if (ownB) { delete B; }
-}
+ProductOperator::ProductOperator(const Operator *A_, const Operator *B_,
+                                 bool own_A, bool own_B)
+   : ProductOperator({A_, own_A}, {B_, own_B}) { }
 
 
 RAPOperator::RAPOperator(const Operator &Rt_, const Operator &A_,
@@ -465,11 +454,9 @@ RAPOperator::RAPOperator(const Operator &Rt_, const Operator &A_,
 
 
 TripleProductOperator::TripleProductOperator(
-   const Operator *A, const Operator *B, const Operator *C,
-   bool ownA, bool ownB, bool ownC)
-   : Operator(A->Height(), C->Width())
-   , A(A), B(B), C(C)
-   , ownA(ownA), ownB(ownB), ownC(ownC)
+   Handle<const Operator> A_, Handle<const Operator> B_, Handle<const Operator> C_)
+   : Operator(A_->Height(), C_->Width()),
+     A(A_), B(B_), C(C_)
 {
    MFEM_VERIFY(A->Width() == B->Height(),
                "incompatible Operators: A->Width() = " << A->Width()
@@ -478,20 +465,16 @@ TripleProductOperator::TripleProductOperator(
                "incompatible Operators: B->Width() = " << B->Width()
                << ", C->Height() = " << C->Height());
 
+   if (auto SolverB = dynamic_cast<const Solver*>(B.Get()))
    {
-      const Solver* SolverB = dynamic_cast<const Solver*>(B);
-      if (SolverB)
-      {
-         MFEM_VERIFY(!(SolverB->iterative_mode),
-                     "Operator B of a TripleProductOperator should not be in iterative mode");
-      }
+      MFEM_VERIFY(!(SolverB->iterative_mode),
+                  "Operator B of a TripleProductOperator should not be in iterative mode");
+   }
 
-      const Solver* SolverC = dynamic_cast<const Solver*>(C);
-      if (SolverC)
-      {
-         MFEM_VERIFY(!(SolverC->iterative_mode),
-                     "Operator C of a TripleProductOperator should not be in iterative mode");
-      }
+   if (auto SolverC = dynamic_cast<const Solver*>(C.Get()))
+   {
+      MFEM_VERIFY(!(SolverC->iterative_mode),
+                  "Operator C of a TripleProductOperator should not be in iterative mode");
    }
 
    mem_class = A->GetMemoryClass()*C->GetMemoryClass();
@@ -500,12 +483,10 @@ TripleProductOperator::TripleProductOperator(
    t2.SetSize(B->Height(), mem_type);
 }
 
-TripleProductOperator::~TripleProductOperator()
-{
-   if (ownA) { delete A; }
-   if (ownB) { delete B; }
-   if (ownC) { delete C; }
-}
+TripleProductOperator::TripleProductOperator(
+   const Operator *A_, const Operator *B_, const Operator *C_,
+   bool own_A, bool own_B, bool own_C)
+   : TripleProductOperator({A_, own_A}, {B_, own_B}, {C_, own_C}) { }
 
 
 ConstrainedOperator::ConstrainedOperator(Operator *A, const Array<int> &list,

linalg/operator.hpp

@@ -13,6 +13,7 @@
 #define MFEM_OPERATOR
 
 #include "vector.hpp"
+#include "../general/handle.hpp"
 
 namespace mfem
 {
@@ -869,43 +870,42 @@ public:
 /// General linear combination operator: x -> a A(x) + b B(x).
 class SumOperator : public Operator
 {
-   const Operator *A, *B;
+   Handle<const Operator> A, B;
    const real_t alpha, beta;
-   bool ownA, ownB;
    mutable Vector z;
 
 public:
-   SumOperator(
-      const Operator *A, const real_t alpha,
-      const Operator *B, const real_t beta,
-      bool ownA, bool ownB);
+   SumOperator(Handle<const Operator> A_, const real_t alpha,
+               Handle<const Operator> B_, const real_t beta);
+
+   SumOperator(const Operator *A_, const real_t alpha,
+               const Operator *B_, const real_t beta,
+               bool own_A, bool own_B);
 
    void Mult(const Vector &x, Vector &y) const override
    { z.SetSize(A->Height()); A->Mult(x, z); B->Mult(x, y); add(alpha, z, beta, y, y); }
 
    void MultTranspose(const Vector &x, Vector &y) const override
    { z.SetSize(A->Width()); A->MultTranspose(x, z); B->MultTranspose(x, y); add(alpha, z, beta, y, y); }
-
-   virtual ~SumOperator();
 };
 
 /// General product operator: x -> (A*B)(x) = A(B(x)).
 class ProductOperator : public Operator
 {
-   const Operator *A, *B;
-   bool ownA, ownB;
+   Handle<const Operator> A, B;
    mutable Vector z;
 
 public:
-   ProductOperator(const Operator *A, const Operator *B, bool ownA, bool ownB);
+   ProductOperator(Handle<const Operator> A_, Handle<const Operator> B_);
+
+   ProductOperator(const Operator *A_, const Operator *B_,
+                   bool own_A, bool own_B);
 
    void Mult(const Vector &x, Vector &y) const override
    { B->Mult(x, z); A->Mult(z, y); }
 
    void MultTranspose(const Vector &x, Vector &y) const override
    { A->MultTranspose(x, z); B->MultTranspose(z, y); }
-
-   virtual ~ProductOperator();
 };
 
 
@@ -956,16 +956,19 @@ public:
 /// General triple product operator x -> A*B*C*x, with ownership of the factors.
 class TripleProductOperator : public Operator
 {
-   const Operator *A;
-   const Operator *B;
-   const Operator *C;
-   bool ownA, ownB, ownC;
+   Handle<const Operator> A;
+   Handle<const Operator> B;
+   Handle<const Operator> C;
    mutable Vector t1, t2;
    MemoryClass mem_class;
 
 public:
-   TripleProductOperator(const Operator *A, const Operator *B,
-                         const Operator *C, bool ownA, bool ownB, bool ownC);
+   TripleProductOperator(Handle<const Operator> A_, Handle<const Operator> B_,
+                         Handle<const Operator> C_);
+
+   TripleProductOperator(
+      const Operator *A_, const Operator *B_, const Operator *C_,
+      bool own_A, bool own_B, bool own_C);
 
    MemoryClass GetMemoryClass() const override { return mem_class; }
 
@@ -974,8 +977,6 @@ public:
 
    void MultTranspose(const Vector &x, Vector &y) const override
    { A->MultTranspose(x, t2); B->MultTranspose(t2, t1); C->MultTranspose(t1, y); }
-
-   virtual ~TripleProductOperator();
 };
 
 

linalg/petsc.hpp

@@ -21,7 +21,7 @@
 
 #include <limits>
 
-#include "handle.hpp"
+#include "op_handle.hpp"
 #include "hypre.hpp"
 #include "ode.hpp"
 #include "../general/mem_manager.hpp"

linalg/solvers.cpp

@@ -1919,7 +1919,8 @@ void NewtonSolver::Mult(const Vector &b, Vector &x) const
        print_options.first_and_last)
    {
       mfem::out << "Newton: Number of iterations: " << final_iter << '\n'
-                << "   ||r|| = " << final_norm << '\n';
+                << "   ||r|| = " << final_norm
+                << ",  ||r||/||r_0|| = " << final_norm/norm0 << '\n';
    }
    if (!converged && (print_options.summary || print_options.warnings))
    {
@@ -2141,7 +2142,8 @@ void LBFGSSolver::Mult(const Vector &b, Vector &x) const
        print_options.first_and_last)
    {
       mfem::out << "LBFGS: Number of iterations: " << final_iter << '\n'
-                << "   ||r|| = " << final_norm << '\n';
+                << "   ||r|| = " << final_norm
+                << ",  ||r||/||r_0|| = " << final_norm/norm0 << '\n';
    }
    if (print_options.summary || (!converged && print_options.warnings))
    {
@@ -2570,7 +2572,7 @@ struct WeightMinHeap
       for (; pos > 0 && w[c[(pos-1)/2]] > val; pos = (pos-1)/2)
       {
          c[pos] = c[(pos-1)/2];
-         loc[c[(pos-1)/2]] = pos;
+         loc[c[(pos-1)/2]] = static_cast<int>(pos);
       }
       return pos;
    }
@@ -2587,7 +2589,7 @@ struct WeightMinHeap
          if (w[c[tgt]] < val)
          {
             c[pos] = c[tgt];
-            loc[c[tgt]] = pos;
+            loc[c[tgt]] = static_cast<int>(pos);
             pos = tgt;
          }
          else
@@ -2605,7 +2607,7 @@ struct WeightMinHeap
       size_t pos = c.size()-1;
       pos = percolate_up(pos, val);
       c[pos] = i;
-      loc[i] = pos;
+      loc[i] = static_cast<int>(pos);
    }
 
    int pop()
@@ -2615,13 +2617,13 @@ struct WeightMinHeap
       c.pop_back();
       // Mark as removed
       loc[i] = -1;
-      if (c.empty()) { return i; }
+      if (c.empty()) { return static_cast<int>(i); }
       real_t val = w[j];
       size_t pos = 0;
       pos = percolate_down(pos, val);
       c[pos] = j;
-      loc[j] = pos;
-      return i;
+      loc[j] = static_cast<int>(pos);
+      return static_cast<int>(i);
    }
 
    void update(size_t i)
@@ -2631,7 +2633,7 @@ struct WeightMinHeap
       pos = percolate_up(pos, val);
       pos = percolate_down(pos, val);
       c[pos] = i;
-      loc[i] = pos;
+      loc[i] = static_cast<int>(pos);
    }
 
    bool picked(size_t i)
@@ -2819,7 +2821,7 @@ void BlockILU::CreateBlockPattern(const SparseMatrix &A)
             unique_block_cols[iblock].insert(J[k] / block_size);
          }
       }
-      nnz += unique_block_cols[iblock].size();
+      nnz += static_cast<int>(unique_block_cols[iblock].size());
    }
 
    if (reordering != Reordering::NONE)

linalg/solvers.hpp

@@ -14,7 +14,7 @@
 
 #include "../config/config.hpp"
 #include "densemat.hpp"
-#include "handle.hpp"
+#include "op_handle.hpp"
 #include <memory>
 
 #ifdef MFEM_USE_MPI

linalg/vector.hpp

@@ -24,6 +24,8 @@
 #include <cstdlib>
 #include <iostream>
 #include <limits>
+#include <type_traits>
+#include <initializer_list>
 #if defined(_MSC_VER) && (_MSC_VER < 1800)
 #include <float.h>
 #define isfinite _finite
@@ -119,10 +121,16 @@ public:
    Vector(int size_, MemoryType h_mt, MemoryType d_mt)
       : data(size_, h_mt, d_mt), size(size_) { }
 
+   /// Create a vector from a statically sized C-style array of convertible type
+   template <typename CT, int N>
+   explicit Vector(const CT (&values)[N]) : Vector(N)
+   { std::copy(values, values + N, begin()); }
+
    /// Create a vector using a braced initializer list
-   template <int N, typename T = real_t>
-   explicit Vector(const T (&values)[N]) : Vector(N)
-   { std::copy(values, values + N, GetData()); }
+   template <typename CT, typename std::enable_if<
+                std::is_convertible<CT,real_t>::value,bool>::type = true>
+   explicit Vector(std::initializer_list<CT> values) : Vector(values.size())
+   { std::copy(values.begin(), values.end(), begin()); }
 
    /// Enable execution of Vector operations using the mfem::Device.
    /** The default is to use Backend::CPU (serial execution on each MPI rank),

mesh/mesh.cpp

@@ -4411,6 +4411,49 @@ Mesh::Mesh(real_t *vertices_, int num_vertices,
    FinalizeTopology();
 }
 
+Mesh::Mesh( const NURBSExtension& ext )
+ : attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
+{
+   SetEmpty();
+   /// make an internal copy of the NURBSExtension
+   NURBSext = new NURBSExtension( ext );
+
+   Dim              = NURBSext->Dimension();
+   NumOfVertices    = NURBSext->GetNV();
+   NumOfElements    = NURBSext->GetNE();
+   NumOfBdrElements = NURBSext->GetNBE();
+
+   NURBSext->GetElementTopo(elements);
+   NURBSext->GetBdrElementTopo(boundary);
+
+   vertices.SetSize(NumOfVertices);
+   if (NURBSext->HavePatches())
+   {
+      NURBSFECollection  *fec = new NURBSFECollection(NURBSext->GetOrder());
+      FiniteElementSpace *fes = new FiniteElementSpace(this, fec, Dim,
+                                                       Ordering::byVDIM);
+      Nodes = new GridFunction(fes);
+      Nodes->MakeOwner(fec);
+      NURBSext->SetCoordsFromPatches(*Nodes);
+      own_nodes = 1;
+      spaceDim = Nodes->VectorDim();
+      for (int i = 0; i < spaceDim; i++)
+      {
+         Vector vert_val;
+         Nodes->GetNodalValues(vert_val, i+1);
+         for (int j = 0; j < NumOfVertices; j++)
+         {
+            vertices[j](i) = vert_val(j);
+         }
+      }
+   }
+   else
+   {
+      MFEM_ABORT("NURBS mesh has no patches.");
+   }
+   FinalizeMesh();
+}
+
 Element *Mesh::NewElement(int geom)
 {
    switch (geom)
@@ -5629,7 +5672,7 @@ std::vector<int> Mesh::CreatePeriodicVertexMapping(
    std::vector<int> v2v(GetNV());
    for (size_t i = 0; i < v2v.size(); i++)
    {
-      v2v[i] = i;
+      v2v[i] = static_cast<int>(i);
    }
    for (const auto &r2p : replica2primary)
    {
@@ -11161,14 +11204,14 @@ const CoarseFineTransformations &Mesh::GetRefinementTransforms() const
             if (code)
             {
                int &matrix = mat_no[code];
-               if (!matrix) { matrix = mat_no.size(); }
+               if (!matrix) { matrix = static_cast<int>(mat_no.size()); }
                index = matrix-1;
             }
             CoarseFineTr.embeddings[j].matrix = index;
          }
 
          DenseTensor &pmats = CoarseFineTr.point_matrices[geom];
-         pmats.SetSize(Dim, Dim+1, mat_no.size());
+         pmats.SetSize(Dim, Dim+1, static_cast<int>((mat_no.size())));
 
          // calculate the point matrices used
          std::map<unsigned, int>::iterator it;
@@ -13888,7 +13931,7 @@ void MeshPartitioner::ExtractPart(int part_id, MeshPart &mesh_part) const
          const int *v = elem->GetVertices();
          vertex_set.insert(v, v + nv);
       }
-      vertex_loc_to_glob.SetSize(vertex_set.size());
+      vertex_loc_to_glob.SetSize(static_cast<int>(vertex_set.size()));
       std::copy(vertex_set.begin(), vertex_set.end(), // src
                 vertex_loc_to_glob.begin());          // dest
    }

mesh/mesh.hpp

@@ -895,6 +895,9 @@ public:
 
    ///@}
 
+   /// Construct a Mesh from a NURBSExtension
+   explicit Mesh( const NURBSExtension& ext );
+
    /** @anchor mfem_Mesh_construction
        @name Methods for piecewise Mesh construction.
 

mesh/mesh_readers.cpp

@@ -746,7 +746,7 @@ struct BufferReader : BufferReaderBase
    /// header_type) pointed to by @a header_buf.
    int NumHeaderBytes(const char *header_buf) const
    {
-      if (!compressed) { return HeaderEntrySize(); }
+      if (!compressed) { return static_cast<int>(HeaderEntrySize()); }
       return (3 + ReadHeaderEntry(header_buf))*HeaderEntrySize();
    }
 
@@ -856,12 +856,13 @@ struct BufferReader : BufferReaderBase
          // Decode the first entry of the header, which we need to determine
          // how long the rest of the header is.
          std::vector<char> nblocks_buf;
-         int nblocks_b64 = bin_io::NumBase64Chars(HeaderEntrySize());
+         int nblocks_b64 = static_cast<int>(bin_io::NumBase64Chars(HeaderEntrySize()));
          bin_io::DecodeBase64(txt, nblocks_b64, nblocks_buf);
          std::vector<char> data, header;
          // Compute number of characters needed to encode header in base 64,
          // then round to nearest multiple of 4 to take padding into account.
-         int header_b64 = bin_io::NumBase64Chars(NumHeaderBytes(nblocks_buf.data()));
+         int header_b64 = static_cast<int>(bin_io::NumBase64Chars(NumHeaderBytes(
+                                                                     nblocks_buf.data())));
          // If data is compressed, header is encoded separately
          bin_io::DecodeBase64(txt, header_b64, header);
          bin_io::DecodeBase64(txt + header_b64, strlen(txt)-header_b64, data);
@@ -2434,13 +2435,13 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
          if (!elements_3D.empty())
          {
             Dim = 3;
-            NumOfElements = elements_3D.size();
+            NumOfElements = static_cast<int>(elements_3D.size());
             elements.SetSize(NumOfElements);
             for (int el = 0; el < NumOfElements; ++el)
             {
                elements[el] = elements_3D[el];
             }
-            NumOfBdrElements = elements_2D.size();
+            NumOfBdrElements = static_cast<int>(elements_2D.size());
             boundary.SetSize(NumOfBdrElements);
             for (int el = 0; el < NumOfBdrElements; ++el)
             {
@@ -2463,13 +2464,13 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
          else if (!elements_2D.empty())
          {
             Dim = 2;
-            NumOfElements = elements_2D.size();
+            NumOfElements = static_cast<int>(elements_2D.size());
             elements.SetSize(NumOfElements);
             for (int el = 0; el < NumOfElements; ++el)
             {
                elements[el] = elements_2D[el];
             }
-            NumOfBdrElements = elements_1D.size();
+            NumOfBdrElements = static_cast<int>(elements_1D.size());
             boundary.SetSize(NumOfBdrElements);
             for (int el = 0; el < NumOfBdrElements; ++el)
             {
@@ -2488,13 +2489,13 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
          else if (!elements_1D.empty())
          {
             Dim = 1;
-            NumOfElements = elements_1D.size();
+            NumOfElements = static_cast<int>(elements_1D.size());
             elements.SetSize(NumOfElements);
             for (int el = 0; el < NumOfElements; ++el)
             {
                elements[el] = elements_1D[el];
             }
-            NumOfBdrElements = elements_0D.size();
+            NumOfBdrElements = static_cast<int>(elements_0D.size());
             boundary.SetSize(NumOfBdrElements);
             for (int el = 0; el < NumOfBdrElements; ++el)
             {

mesh/ncmesh.cpp

@@ -3084,13 +3084,13 @@ struct MatrixMap
    int GetIndex(const NCMesh::PointMatrix &pm)
    {
       int &index = map[pm];
-      if (!index) { index = map.size(); }
+      if (!index) { index = static_cast<int>(map.size()); }
       return index - 1;
    }
 
    void ExportMatrices(Array<DenseMatrix*> &point_matrices) const
    {
-      point_matrices.SetSize(map.size());
+      point_matrices.SetSize(static_cast<int>(map.size()));
       for (const auto &pair : map)
       {
          DenseMatrix* mat = new DenseMatrix();
@@ -4813,7 +4813,7 @@ void NCMesh::TraverseRefinements(int elem, int coarse_index,
    if (!el.ref_type)
    {
       int &matrix = map[ref_path];
-      if (!matrix) { matrix = map.size(); }
+      if (!matrix) { matrix = static_cast<int>(map.size()); }
 
       Embedding &emb = transforms.embeddings[el.index];
       emb.parent = coarse_index;
@@ -4876,7 +4876,7 @@ const CoarseFineTransformations& NCMesh::GetRefinementTransforms() const
             const PointMatrix &identity = GetGeomIdentity(geom);
 
             transforms.point_matrices[g]
-            .SetSize(Dim, identity.np, path_map[g].size());
+            .SetSize(Dim, identity.np, static_cast<int>(path_map[g].size()));
 
             // calculate the point matrices
             RefPathMap::iterator it;
@@ -4912,7 +4912,7 @@ const CoarseFineTransformations& NCMesh::GetDerefinementTransforms() const
          if (code)
          {
             int &matrix = mat_no[emb.geom][code];
-            if (!matrix) { matrix = mat_no[emb.geom].size(); }
+            if (!matrix) { matrix = static_cast<int>(mat_no[emb.geom].size()); }
 
             emb.matrix = matrix - 1;
          }
@@ -4926,7 +4926,7 @@ const CoarseFineTransformations& NCMesh::GetDerefinementTransforms() const
             const PointMatrix &identity = GetGeomIdentity(geom);
 
             transforms.point_matrices[geom]
-            .SetSize(Dim, identity.np, mat_no[geom].size());
+            .SetSize(Dim, identity.np, static_cast<int>(mat_no[geom].size()));
 
             // calculate point matrices
             for (auto it = mat_no[geom].begin(); it != mat_no[geom].end(); ++it)
@@ -6529,9 +6529,9 @@ long NCMesh::NCList::MemoryUsage() const
    {
       for (int j = 0; j < point_matrices[i].Size(); i++)
       {
-         pm_size += point_matrices[i][j]->MemoryUsage();
+         pm_size += static_cast<int>(point_matrices[i][j]->MemoryUsage());
       }
-      pm_size += point_matrices[i].MemoryUsage();
+      pm_size += static_cast<int>(point_matrices[i].MemoryUsage());
    }
 
    return conforming.MemoryUsage() +

mesh/nurbs.cpp

@@ -45,6 +45,37 @@ KnotVector::KnotVector(int order, int NCP)
    knot = -1.;
 }
 
+KnotVector::KnotVector(int order, const Vector& intervals,
+                       const Array<int>& continuity )
+{
+   // NOTE: This may need to be generalized to support periodicity
+   // in the future.
+   MFEM_ASSERT(continuity.Size() == (intervals.Size() + 1),
+               "Incompatible sizes of continuity and intervals.");
+   Order = order;
+   NumOfElements = intervals.Size();
+   const int num_knots = Order * continuity.Size() - continuity.Sum();
+   // Some continuities may still be invalid; this assert only avoids
+   // passing a negative num_knots to Vector::SetSize().
+   MFEM_ASSERT(num_knots >= 0, "Invalid continuity vector for order.");
+   NumOfControlPoints = num_knots - Order - 1;
+   knot.SetSize(num_knots);
+   real_t accum = 0.0;
+   int iknot = 0;
+   for (int i = 0; i < continuity.Size(); ++i)
+   {
+      const int multiplicity = Order - continuity[i];
+      MFEM_ASSERT(multiplicity >= 1 && multiplicity <= Order + 1,
+                  "Invalid knot multiplicity for order.");
+      for (int j = 0; j < multiplicity; ++j)
+      {
+         knot[iknot] = accum;
+         ++iknot;
+      }
+      if (i < intervals.Size()) { accum += intervals[i]; }
+   }
+}
+
 KnotVector &KnotVector::operator=(const KnotVector &kv)
 {
    Order = kv.Order;
@@ -54,9 +85,6 @@ KnotVector &KnotVector::operator=(const KnotVector &kv)
    coarse = kv.coarse;
    if (kv.spacing) { spacing = kv.spacing->Clone(); }
 
-   // alternatively, re-compute NumOfElements
-   // GetElements();
-
    return *this;
 }
 
@@ -640,12 +668,15 @@ void KnotVector::Difference(const KnotVector &kv, Vector &diff) const
 
 void NURBSPatch::init(int dim)
 {
+   MFEM_ASSERT(dim > 1, "NURBS patch dimension (including weight) must be "
+               "greater than 1.");
    Dim = dim;
    sd = nd = -1;
 
    if (kv.Size() == 1)
    {
       ni = kv[0]->GetNCP();
+      MFEM_ASSERT(ni > 0, "Invalid knot vector dimension.");
       nj = -1;
       nk = -1;
 
@@ -662,6 +693,7 @@ void NURBSPatch::init(int dim)
    {
       ni = kv[0]->GetNCP();
       nj = kv[1]->GetNCP();
+      MFEM_ASSERT(ni > 0 && nj > 0, "Invalid knot vector dimensions.");
       nk = -1;
 
       data = new real_t[ni*nj*Dim];
@@ -678,6 +710,8 @@ void NURBSPatch::init(int dim)
       ni = kv[0]->GetNCP();
       nj = kv[1]->GetNCP();
       nk = kv[2]->GetNCP();
+      MFEM_ASSERT(ni > 0 && nj > 0 && nk > 0,
+                  "Invalid knot vector dimensions.");
 
       data = new real_t[ni*nj*nk*Dim];
 
@@ -2330,13 +2364,77 @@ NURBSExtension::NURBSExtension(Mesh *mesh_array[], int num_pieces)
    MergeWeights(mesh_array, num_pieces);
 }
 
+NURBSExtension::NURBSExtension(const Mesh *patch_topology,
+                               const Array<const NURBSPatch*> p)
+{
+   patchTopo = new Mesh( *patch_topology );
+   patchTopo->GetEdgeVertexTable();
+   own_topo = 1;
+   patches.Reserve(p.Size());
+   Array<int> edges;
+   Array<int> oedges;
+   Array<int> kvs(3);
+   edge_to_knot.SetSize(patch_topology->GetNEdges());
+   NumOfKnotVectors = 0;
+   NumOfElements = 0;
+   for (int ielem = 0; ielem < patch_topology->GetNE(); ++ielem)
+   {
+      patches.Append(new NURBSPatch(*p[ielem]));
+      NURBSPatch& patch = *patches[ielem];
+      int num_patch_elems = 1;
+      for (int ikv = 0; ikv < patch.GetNKV(); ++ikv)
+      {
+         kvs[ikv] = knotVectors.Size();
+         knotVectors.Append(new KnotVector(*patch.GetKV(ikv)));
+         num_patch_elems *= patch.GetKV(ikv)->GetNE();
+         ++NumOfKnotVectors;
+      }
+      NumOfElements += num_patch_elems;
+      patch_topology->GetElementEdges(ielem, edges, oedges);
+      for (int iedge = 0; iedge < edges.Size(); ++iedge)
+      {
+         if (iedge < 8)
+         {
+            if (iedge & 1)
+            {
+               edge_to_knot[edges[iedge]] = kvs[1];
+            }
+            else
+            {
+               edge_to_knot[edges[iedge]] = kvs[0];
+            }
+         }
+         else
+         {
+            edge_to_knot[edges[iedge]] = kvs[2];
+         }
+      }
+   }
+
+   GenerateOffsets();
+   CountBdrElements();
+   NumOfActiveElems = NumOfElements;
+   activeElem.SetSize(NumOfElements);
+   activeElem = true;
+
+   CreateComprehensiveKV();
+   SetOrdersFromKnotVectors();
+
+   GenerateActiveVertices();
+   InitDofMap();
+   GenerateElementDofTable();
+   GenerateActiveBdrElems();
+   GenerateBdrElementDofTable();
+
+   weights.SetSize(GetNDof());
+
+   CheckPatches();
+}
+
 NURBSExtension::~NURBSExtension()
 {
-   if (patches.Size() == 0)
-   {
-      delete bel_dof;
-      delete el_dof;
-   }
+   if (bel_dof) { delete bel_dof; }
+   if (el_dof) { delete el_dof; }
 
    for (int i = 0; i < knotVectors.Size(); i++)
    {
@@ -4804,6 +4902,42 @@ const Array<int>& NURBSExtension::GetPatchBdrElements(int patch)
    return patch_to_bel[patch];
 }
 
+NURBSPatch::NURBSPatch(const KnotVector *kv0, const KnotVector *kv1, int dim_,
+                       const real_t* control_points)
+{
+   kv.SetSize(2);
+   kv[0] = new KnotVector(*kv0);
+   kv[1] = new KnotVector(*kv1);
+   init(dim_);
+   memcpy(data, control_points, sizeof (real_t) * ni * nj * dim_);
+}
+
+NURBSPatch::NURBSPatch(const KnotVector *kv0, const KnotVector *kv1,
+                       const KnotVector *kv2, int dim_,
+                       const real_t* control_points)
+{
+   kv.SetSize(3);
+   kv[0] = new KnotVector(*kv0);
+   kv[1] = new KnotVector(*kv1);
+   kv[2] = new KnotVector(*kv2);
+   init(dim_);
+   memcpy(data, control_points, sizeof (real_t) * ni * nj * nk * dim_);
+}
+
+NURBSPatch::NURBSPatch(Array<const KnotVector *> &kv_,  int dim_,
+                       const real_t* control_points)
+{
+   kv.SetSize(kv_.Size());
+   int n = dim_;
+   for (int i = 0; i < kv.Size(); i++)
+   {
+      kv[i] = new KnotVector(*kv_[i]);
+      n *= kv[i]->GetNCP();
+   }
+   init(dim_);
+   memcpy(data, control_points, sizeof(real_t)*n);
+}
+
 #ifdef MFEM_USE_MPI
 ParNURBSExtension::ParNURBSExtension(const ParNURBSExtension &orig)
    : NURBSExtension(orig),

mesh/nurbs.hpp

@@ -65,6 +65,17 @@ public:
        order @a order and number of control points @a NCP. */
    KnotVector(int order, int NCP);
 
+   /** @brief Create a KnotVector by passing in a degree, a Vector of interval
+       lengths of length n, and a list of continuity of length n + 1.
+
+       The intervals refer to spans between unique knot values (not counting
+       zero-size intervals at repeated knots), and the continuity values should
+       be >= -1 (discontinuous) and <= order-1 (maximally-smooth for the given
+       polynomial degree). Periodicity is not supported.
+   */
+   KnotVector(int order, const Vector& intervals,
+              const Array<int>& continuity );
+
    /// Copy constructor.
    KnotVector(const KnotVector &kv) { (*this) = kv; }
 
@@ -262,6 +273,22 @@ public:
    NURBSPatch(const KnotVector *kv0, const KnotVector *kv1,
               const KnotVector *kv2, int dim);
 
+   /** Create a bivariate NURBS patch with given control points. See n-variate
+       overload for additional notes. */
+   NURBSPatch(const KnotVector *kv0, const KnotVector *kv1, int dim_,
+              const real_t* control_points);
+   /** Create a trivariate NURBS patch with given control points. See n-variate
+       overload for additional notes. */
+   NURBSPatch(const KnotVector *kv0, const KnotVector *kv1,
+              const KnotVector *kv2, int dim_, const real_t* control_points);
+   /** Create an n-variate NURBS patch with given control points of dimension
+       dim_, where n is the length of the array of knot vectors and dim_
+       includes the weight. The array of control point coordinates stores each
+       point's coordinates contiguously, and points are ordered in a standard
+       ijk grid ordering. */
+   NURBSPatch(Array<const KnotVector *> &kv_,  int dim_,
+              const real_t* control_points);
+
    /// Constructor for a patch of dimension equal to the size of @a kv.
    NURBSPatch(Array<const KnotVector *> &kv, int dim);
 
@@ -647,7 +674,7 @@ protected:
    void SetPatchToBdrElements();
 
    /// To be used by ParNURBSExtension constructor(s)
-   NURBSExtension() { }
+   NURBSExtension() : el_dof(nullptr), bel_dof(nullptr) { }
 
 public:
    /// Copy constructor: deep copy
@@ -671,6 +698,8 @@ public:
 
    NURBSExtension(Mesh *mesh_array[], int num_pieces);
 
+   NURBSExtension(const Mesh *patch_topology, const Array<const NURBSPatch*> p);
+
    /// Copy assignment not supported.
    NURBSExtension& operator=(const NURBSExtension&) = delete;
 

mesh/pncmesh.cpp

@@ -741,7 +741,7 @@ void ParNCMesh::ElementNeighborProcessors(int elem, Array<int> &ranks)
 template<class T>
 static void set_to_array(const std::set<T> &set, Array<T> &array)
 {
-   array.Reserve(set.size());
+   array.Reserve(static_cast<int>(set.size()));
    array.SetSize(0);
    for (auto x : set)
    {
@@ -842,7 +842,7 @@ void ParNCMesh::GetConformingSharedStructures(ParMesh &pmesh)
    }
 
    // create ParMesh groups, and the map (ncmesh_group -> pmesh_group)
-   Array<int> group_map(groups.size());
+   Array<int> group_map(static_cast<int>(groups.size()));
    {
       group_map = 0;
       IntegerSet iset;
@@ -851,7 +851,7 @@ void ParNCMesh::GetConformingSharedStructures(ParMesh &pmesh)
       {
          if (groups[i].size() > 1 || !i) // skip singleton groups
          {
-            iset.Recreate(groups[i].size(), groups[i].data());
+            iset.Recreate(static_cast<int>(groups[i].size()), groups[i].data());
             group_map[i] = int_groups.Insert(iset);
          }
       }
@@ -1052,7 +1052,7 @@ void ParNCMesh::GetFaceNeighbors(ParMesh &pmesh)
       for (int k = 0; k < gi.nv; k++)
       {
          int &v = vert_map[elem->node[k]];
-         if (!v) { v = vert_map.size(); }
+         if (!v) { v = static_cast<int>(vert_map.size()); }
          fne->GetVertices()[k] = v-1;
       }
 
@@ -1068,7 +1068,7 @@ void ParNCMesh::GetFaceNeighbors(ParMesh &pmesh)
 
    // create vertices in 'face_nbr_vertices'
    {
-      pmesh.face_nbr_vertices.SetSize(vert_map.size());
+      pmesh.face_nbr_vertices.SetSize(static_cast<int>(vert_map.size()));
       if (coordinates.Size())
       {
          tmp_vertex = new TmpVertex[nodes.NumIds()]; // TODO: something cheaper?
@@ -2231,7 +2231,7 @@ void ParNCMesh::SendRebalanceDofs(int old_ndofs,
    {
       RebalanceDofMessage &msg = it->second;
       msg.dofs.clear();
-      int ne = msg.elem_ids.size();
+      int ne = static_cast<int>(msg.elem_ids.size());
       if (ne)
       {
          msg.dofs.reserve(old_element_dofs.RowSize(msg.elem_ids[0]) * ne * vdim);
@@ -2261,8 +2261,8 @@ void ParNCMesh::RecvRebalanceDofs(Array<int> &elements, Array<long> &dofs)
    for (it = recv_rebalance_dofs.begin(); it != recv_rebalance_dofs.end(); ++it)
    {
       RebalanceDofMessage &msg = it->second;
-      ne += msg.elem_ids.size();
-      nd += msg.dofs.size();
+      ne += static_cast<int>(msg.elem_ids.size());
+      nd += static_cast<int>(msg.dofs.size());
    }
 
    elements.SetSize(ne);
@@ -2480,7 +2480,7 @@ void ParNCMesh::AdjustMeshIds(Array<MeshId> ids[], int rank)
    if (!shared_edges.masters.Size() &&
        !shared_faces.masters.Size()) { return; }
 
-   Array<bool> contains_rank(groups.size());
+   Array<bool> contains_rank(static_cast<int>(groups.size()));
    for (unsigned i = 0; i < groups.size(); i++)
    {
       contains_rank[i] = GroupContains(i, rank);
@@ -2795,7 +2795,7 @@ void ParNCMesh::ElementValueMessage<ValueType, RefTypes, Tag>::Encode(int)
    std::ostringstream ostream;
 
    Array<int> tmp_elements;
-   tmp_elements.MakeRef(elements.data(), elements.size());
+   tmp_elements.MakeRef(elements.data(), static_cast<int>(elements.size()));
 
    ElementSet eset(pncmesh, RefTypes);
    eset.Encode(tmp_elements);
@@ -2812,7 +2812,7 @@ void ParNCMesh::ElementValueMessage<ValueType, RefTypes, Tag>::Encode(int)
       element_index[decoded[i]] = i;
    }
 
-   write<int>(ostream, values.size());
+   write<int>(ostream, static_cast<int>(values.size()));
    MFEM_ASSERT(elements.size() == values.size(), "");
 
    for (unsigned i = 0; i < values.size(); i++)
@@ -2870,7 +2870,7 @@ void ParNCMesh::RebalanceDofMessage::SetElements(const Array<int> &elems,
 
 static void write_dofs(std::ostream &os, const std::vector<int> &dofs)
 {
-   write<int>(os, dofs.size());
+   write<int>(os, static_cast<int>(dofs.size()));
    // TODO: we should compress the ints, mostly they are contiguous ranges
    os.write((const char*) dofs.data(), dofs.size() * sizeof(int));
 }

mesh/pncmesh.hpp

@@ -441,7 +441,7 @@ protected: // implementation
       std::vector<int> elements;
       std::vector<ValueType> values;
 
-      int Size() const { return elements.size(); }
+      int Size() const { return static_cast<int>(elements.size()); }
       void Reserve(int size) { elements.reserve(size); values.reserve(size); }
 
       void Add(int elem, ValueType val)

mesh/submesh/submesh_utils.cpp

@@ -611,7 +611,8 @@ void ConstructFaceTree(NCSubMeshT &submesh, const Array<int> &attributes)
          stable across processors.
    */
    // Build an inverse (and consecutive) map.
-   Array<FaceNodes> new_elem_to_parent_face_nodes(pnodes_new_elem.size());
+   Array<FaceNodes> new_elem_to_parent_face_nodes(static_cast<int>
+                                                  (pnodes_new_elem.size()));
    for (const auto &kv : pnodes_new_elem)
    {
       new_elem_to_parent_face_nodes[kv.second] = kv.first;

miniapps/dpg/convection-diffusion.cpp

@@ -68,7 +68,7 @@
 #include <fstream>
 #include <iostream>
 
-using namespace std;
+
 using namespace mfem;
 using namespace mfem::common;
 
@@ -131,7 +131,7 @@ int main(int argc, char *argv[])
    args.Parse();
    if (!args.Good())
    {
-      args.PrintUsage(cout);
+      args.PrintUsage(std::cout);
       return 1;
    }
 
@@ -154,7 +154,7 @@ int main(int argc, char *argv[])
       beta[0] = 1.;
    }
 
-   args.PrintOptions(cout);
+   args.PrintOptions(std::cout);
 
    // Define spaces
    enum TrialSpace
@@ -292,8 +292,8 @@ int main(int argc, char *argv[])
              << "  L2 Error  |"
              << "  Rate  |"
              << "  Residual  |"
-             << "  Rate  |" << endl;
-   std::cout << std::string(64,'-') << endl;
+             << "  Rate  |" << std::endl;
+   std::cout << std::string(64,'-') << std::endl;
 
    if (static_cond) { a->EnableStaticCondensation(); }
    for (int it = 0; it<=ref; it++)
@@ -625,8 +625,8 @@ void setup_test_norm_coeffs(GridFunction & c1_gf, GridFunction & c2_gf)
    for (int i = 0; i < mesh->GetNE(); i++)
    {
       real_t volume = mesh->GetElementVolume(i);
-      real_t c1 = min(epsilon/volume, (real_t) 1.);
-      real_t c2 = min(1./epsilon, 1./volume);
+      real_t c1 = std::min(epsilon/volume, (real_t) 1.);
+      real_t c2 = std::min(1./epsilon, 1./volume);
       fes->GetElementDofs(i,vdofs);
       c1_gf.SetSubVector(vdofs,c1);
       c2_gf.SetSubVector(vdofs,c2);

miniapps/dpg/pconvection-diffusion.cpp

@@ -70,7 +70,7 @@
 #include <fstream>
 #include <iostream>
 
-using namespace std;
+
 using namespace mfem;
 using namespace mfem::common;
 
@@ -157,7 +157,7 @@ int main(int argc, char *argv[])
    {
       if (myid == 0)
       {
-         args.PrintUsage(cout);
+         args.PrintUsage(std::cout);
       }
       return 1;
    }
@@ -204,7 +204,7 @@ int main(int argc, char *argv[])
 
    if (myid == 0)
    {
-      args.PrintOptions(cout);
+      args.PrintOptions(std::cout);
    }
 
    mesh.EnsureNCMesh(true);
@@ -353,14 +353,14 @@ int main(int argc, char *argv[])
                 << "    Dofs    |" ;
       if (exact_known)
       {
-         mfem::out << "  L2 Error  |"
+         std::cout << "  L2 Error  |"
                    << "  Rate  |";
       }
       std::cout << "  Residual  |"
                 << "  Rate  |"
-                << " CG it  |" << endl;
+                << " CG it  |" << std::endl;
       std::cout << std::string((exact_known) ? 72 : 50,'-')
-                << endl;
+                << std::endl;
    }
 
    if (static_cond) { a->EnableStaticCondensation(); }
@@ -858,8 +858,8 @@ void setup_test_norm_coeffs(ParGridFunction & c1_gf, ParGridFunction & c2_gf)
    for (int i = 0; i < pmesh->GetNE(); i++)
    {
       real_t volume = pmesh->GetElementVolume(i);
-      real_t c1 = min(epsilon/volume, (real_t) 1.);
-      real_t c2 = min(1./epsilon, 1./volume);
+      real_t c1 = std::min(epsilon/volume, (real_t) 1.);
+      real_t c2 = std::min(1./epsilon, 1./volume);
       fes->GetElementDofs(i,vdofs);
       c1_gf.SetSubVector(vdofs,c1);
       c2_gf.SetSubVector(vdofs,c2);

miniapps/meshing/mesh-optimizer.cpp

@@ -51,14 +51,13 @@
 //   Adapted discrete size NC mesh;
 //     mesh-optimizer -m amr-quad-q2.mesh -o 2 -rs 2 -mid 94 -tid 5 -ni 50 -qo 4 -nor
 //   Adapted discrete size 3D with PA:
-//     mesh-optimizer -m cube.mesh -o 2 -rs 2 -mid 321 -tid 5 -ls 3 -nor -pa
+//     mesh-optimizer -m cube.mesh -o 2 -rs 2 -mid 321 -tid 5 -ls 3 -nor -pa -rtol 1e-8
 //   Adapted discrete size 3D with PA on device (requires CUDA):
 //   * mesh-optimizer -m cube.mesh -o 3 -rs 3 -mid 321 -tid 5 -ls 3 -nor -lc 0.1 -pa -d cuda
 //   Adapted discrete size; explicit combo of metrics; mixed tri/quad mesh:
 //     mesh-optimizer -m ../../data/square-mixed.mesh -o 2 -rs 2 -mid 2 -tid 5 -ni 200 -bnd -qo 6 -cmb 2 -nor
 //   Adapted discrete size+aspect_ratio:
 //     mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 6 -ni 100
-//     mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 6 -ni 100 -qo 6 -ex -st 1 -nor
 //   Adapted discrete size+orientation:
 //      mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 36 -tid 8 -qo 4 -fd -nor
 //   Adapted discrete aspect ratio (3D):
@@ -67,7 +66,7 @@
 //   Adaptive limiting:
 //     mesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 50 -qo 5 -nor -vl 1 -alc 0.5
 //   Adaptive limiting through the L-BFGS solver:
-//     mesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 400 -qo 5 -nor -vl 1 -alc 0.5 -st 1
+//     mesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 400 -qo 5 -nor -vl 1 -alc 0.5 -st 1 -rtol 1e-8
 //   Adaptive limiting through FD (requires GSLIB):
 //   * mesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 50 -qo 5 -nor -vl 1 -alc 0.5 -fd -ae 1
 //
@@ -558,7 +557,9 @@ int main(int argc, char *argv[])
    TargetConstructor *target_c = NULL;
    HessianCoefficient *adapt_coeff = NULL;
    HRHessianCoefficient *hr_adapt_coeff = NULL;
-   H1_FECollection ind_fec(mesh_poly_deg, dim);
+   int ind_fec_order = (target_id >= 5 && target_id <= 8 && !fdscheme) ?
+                       1 : mesh_poly_deg;
+   H1_FECollection ind_fec(ind_fec_order, dim);
    FiniteElementSpace ind_fes(mesh, &ind_fec);
    FiniteElementSpace ind_fesv(mesh, &ind_fec, dim);
    GridFunction size(&ind_fes), aspr(&ind_fes), ori(&ind_fes);
@@ -599,6 +600,7 @@ int main(int argc, char *argv[])
          }
          ConstructSizeGF(size);
          tc->SetSerialDiscreteTargetSize(size);
+         tc->SetMinSizeForTargets(size.Min());
          target_c = tc;
          break;
       }
@@ -693,6 +695,7 @@ int main(int argc, char *argv[])
          DiffuseField(aspr, 2);
 
          tc->SetSerialDiscreteTargetSize(size);
+         tc->SetMinSizeForTargets(size.Min());
          tc->SetSerialDiscreteTargetAspectRatio(aspr);
          target_c = tc;
          break;
@@ -740,6 +743,7 @@ int main(int argc, char *argv[])
          ConstantCoefficient size_coeff(0.1*0.1);
          size.ProjectCoefficient(size_coeff);
          tc->SetSerialDiscreteTargetSize(size);
+         tc->SetMinSizeForTargets(size.Min());
 
          FunctionCoefficient ori_coeff(discrete_ori_2d);
          ori.ProjectCoefficient(ori_coeff);
@@ -1109,8 +1113,8 @@ int main(int argc, char *argv[])
    }
    // Level of output.
    IterativeSolver::PrintLevel newton_print;
-   if (verbosity_level > 0)
-   { newton_print.Errors().Warnings().Iterations(); }
+   if (verbosity_level > 0) { newton_print.Errors().Warnings().Iterations(); }
+   else { newton_print.Errors().Warnings(); }
    solver.SetPrintLevel(newton_print);
    // hr-adaptivity solver.
    // If hr-adaptivity is disabled, r-adaptivity is done once using the

miniapps/meshing/pmesh-fitting.cpp

@@ -37,11 +37,11 @@
 //    mpirun -np 4 pmesh-fitting -o 3 -mid 58 -tid 1 -vl 1 -sfc 5e4 -rtol 1e-5
 //    mpirun -np 4 pmesh-fitting -m square01-tri.mesh -o 3 -rs 0 -mid 58 -tid 1 -vl 1 -sfc 1e4 -rtol 1e-5
 //  Surface fitting with weight adaptation and termination based on fitting error:
-//    mpirun -np 4 pmesh-fitting -o 2 -mid 2 -tid 1 -vl 2 -sfc 10 -rtol 1e-20 -sfa 10.0 -sft 1e-5 -no-resid
+//    mpirun -np 4 pmesh-fitting -o 2 -mid 2 -tid 1 -vl 2 -sfc 10 -rtol 1e-20 -sfa 10.0 -sft 1e-5 -no-resid -ni 40
 //  Surface fitting with weight adaptation, limit on max weight, and convergence based on residual.
 //  * mpirun -np 4 pmesh-fitting -m ../../data/inline-tri.mesh -o 2 -mid 2 -tid 4 -vl 2 -sfc 10 -rtol 1e-10 -sfa 10.0 -sft 1e-5 -bgamriter 3 -sbgmesh -ae 1 -marking -slstype 3 -resid -sfcmax 10000 -mod-bndr-attr
 //  Surface fitting to Fischer-Tropsch reactor like domain (requires GSLIB):
-//  * mpirun -np 6 pmesh-fitting -m ../../data/inline-tri.mesh -o 2 -rs 4 -mid 2 -tid 1 -vl 2 -sfc 100 -rtol 1e-12 -li 20 -ae 1 -bnd -sbgmesh -slstype 2 -smtype 0 -sfa 10.0 -sft 1e-4 -no-resid -bgamriter 5 -dist -mod-bndr-attr
+//  * mpirun -np 6 pmesh-fitting -m ../../data/inline-tri.mesh -o 2 -rs 4 -mid 2 -tid 1 -vl 2 -sfc 100 -rtol 1e-12 -li 20 -ae 1 -bnd -sbgmesh -slstype 2 -smtype 0 -sfa 10.0 -sft 1e-4 -no-resid -bgamriter 5 -dist -mod-bndr-attr -ni 50
 
 #include "mesh-fitting.hpp"
 
@@ -809,7 +809,7 @@ int main (int argc, char *argv[])
    solver.SetRelTol(solver_rtol);
    solver.SetAbsTol(0.0);
    solver.SetMinimumDeterminantThreshold(0.001*min_detJ);
-   solver.SetPrintLevel(verbosity_level >= 1 ? 1 : -1);
+   solver.SetPrintLevel(verbosity_level >= 1 ? 1 : 0);
    solver.SetOperator(a);
    Vector b(0);
    solver.Mult(b, x.GetTrueVector());

miniapps/meshing/pmesh-optimizer.cpp

@@ -51,14 +51,13 @@
 //   Adapted discrete size NC mesh;
 //     mpirun -np 4 pmesh-optimizer -m amr-quad-q2.mesh -o 2 -rs 2 -mid 94 -tid 5 -ni 50 -qo 4 -nor
 //   Adapted discrete size 3D with PA:
-//     mpirun -np 4 pmesh-optimizer -m cube.mesh -o 2 -rs 2 -mid 321 -tid 5 -ls 3 -nor -pa
+//     mpirun -np 4 pmesh-optimizer -m cube.mesh -o 2 -rs 2 -mid 321 -tid 5 -ls 3 -nor -pa -rtol 1e-8
 //   Adapted discrete size 3D with PA on device (requires CUDA):
 //   * mpirun -n 4 pmesh-optimizer -m cube.mesh -o 3 -rs 3 -mid 321 -tid 5 -ls 3 -nor -lc 0.1 -pa -d cuda
 //   Adapted discrete size; explicit combo of metrics; mixed tri/quad mesh:
 //     mpirun -np 4 pmesh-optimizer -m ../../data/square-mixed.mesh -o 2 -rs 2 -mid 2 -tid 5 -ni 200 -bnd -qo 6 -cmb 2 -nor
 //   Adapted discrete size+aspect_ratio:
 //     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 6 -ni 100
-//     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 6 -ni 100 -qo 6 -ex -st 1 -nor
 //   Adapted discrete size+orientation:
 //     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 36 -tid 8 -qo 4 -fd -nor
 //   Adapted discrete aspect ratio (3D):
@@ -67,7 +66,7 @@
 //   Adaptive limiting:
 //     mpirun -np 4 pmesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 50 -qo 5 -nor -vl 1 -alc 0.5
 //   Adaptive limiting through the L-BFGS solver:
-//     mpirun -np 4 pmesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 400 -qo 5 -nor -vl 1 -alc 0.5 -st 1
+//     mpirun -np 4 pmesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 400 -qo 5 -nor -vl 1 -alc 0.5 -st 1 -rtol 1e-8
 //   Adaptive limiting through FD (requires GSLIB):
 //   * mpirun -np 4 pmesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 50 -qo 5 -nor -vl 1 -alc 0.5 -fd -ae 1
 //
@@ -583,7 +582,9 @@ int main (int argc, char *argv[])
    TargetConstructor *target_c = NULL;
    HessianCoefficient *adapt_coeff = NULL;
    HRHessianCoefficient *hr_adapt_coeff = NULL;
-   H1_FECollection ind_fec(mesh_poly_deg, dim);
+   int ind_fec_order = (target_id >= 5 && target_id <= 8 && !fdscheme) ?
+                       1 : mesh_poly_deg;
+   H1_FECollection ind_fec(ind_fec_order, dim);
    ParFiniteElementSpace ind_fes(pmesh, &ind_fec);
    ParFiniteElementSpace ind_fesv(pmesh, &ind_fec, dim);
    ParGridFunction size(&ind_fes), aspr(&ind_fes), ori(&ind_fes);
@@ -624,6 +625,7 @@ int main (int argc, char *argv[])
          }
          ConstructSizeGF(size);
          tc->SetParDiscreteTargetSize(size);
+         tc->SetMinSizeForTargets(size.Min());
          target_c = tc;
          break;
       }
@@ -726,6 +728,7 @@ int main (int argc, char *argv[])
          DiffuseField(aspr, 2);
 
          tc->SetParDiscreteTargetSize(size);
+         tc->SetMinSizeForTargets(size.Min());
          tc->SetParDiscreteTargetAspectRatio(aspr);
          target_c = tc;
          break;
@@ -772,6 +775,7 @@ int main (int argc, char *argv[])
          ConstantCoefficient size_coeff(0.1*0.1);
          size.ProjectCoefficient(size_coeff);
          tc->SetParDiscreteTargetSize(size);
+         tc->SetMinSizeForTargets(size.Min());
 
          FunctionCoefficient ori_coeff(discrete_ori_2d);
          ori.ProjectCoefficient(ori_coeff);
@@ -1154,8 +1158,8 @@ int main (int argc, char *argv[])
    }
    // Level of output.
    IterativeSolver::PrintLevel newton_print;
-   if (verbosity_level > 0)
-   { newton_print.Errors().Warnings().Iterations(); }
+   if (verbosity_level > 0) { newton_print.Errors().Warnings().Iterations(); }
+   else { newton_print.Errors().Warnings(); }
    solver.SetPrintLevel(newton_print);
    // hr-adaptivity solver.
    // If hr-adaptivity is disabled, r-adaptivity is done once using the

tests/unit/fem/test_pa_kernels.cpp

@@ -830,3 +830,33 @@ TEST_CASE("Parallel PA DG Diffusion", "[PartialAssembly][Parallel][CUDA]")
 #endif
 
 } // namespace pa_kernels
+
+TEST_CASE("Dispatch Map Specializations")
+{
+   // The kernel specializations are registered the first time the associated
+   // object is created (in the constructor of a static local variable in the
+   // object's constructor). We create a dummy objects here to ensure that the
+   // kernels are registered before testing.
+
+   MassIntegrator{};
+   REQUIRE_FALSE(MassIntegrator::ApplyPAKernels::GetDispatchTable().empty());
+   REQUIRE_FALSE(MassIntegrator::DiagonalPAKernels::GetDispatchTable().empty());
+
+   DiffusionIntegrator{};
+   REQUIRE_FALSE(
+      DiffusionIntegrator::ApplyPAKernels::GetDispatchTable().empty());
+   REQUIRE_FALSE(
+      DiffusionIntegrator::DiagonalPAKernels::GetDispatchTable().empty());
+
+   Mesh mesh = Mesh::MakeCartesian2D(2, 2, Element::QUADRILATERAL);
+   H1_FECollection fec(1, mesh.Dimension());
+   FiniteElementSpace fes(&mesh, &fec);
+   fes.GetQuadratureInterpolator(IntRules.Get(mesh.GetElementGeometry(0), 1));
+
+   using QI = QuadratureInterpolator;
+   REQUIRE_FALSE(QI::TensorEvalKernels::GetDispatchTable().empty());
+   REQUIRE_FALSE(QI::GradKernels::GetDispatchTable().empty());
+   REQUIRE_FALSE(QI::DetKernels::GetDispatchTable().empty());
+   REQUIRE_FALSE(QI::EvalKernels::GetDispatchTable().empty());
+   REQUIRE_FALSE(QI::CollocatedGradKernels::GetDispatchTable().empty());
+}

tests/unit/fem/test_pgridfunc_save_serial.cpp

@@ -14,79 +14,73 @@
 
 using namespace mfem;
 
-namespace pgridfunc_save_in_serial
-{
-
-double squared(const Vector &x)
-{
-   double sum = 0.0;
-   for (int d = 0; d < x.Size(); d++)
-   {
-      sum += std::pow(x(d), 2.0);
-   }
-   return sum;
-}
-
-
 #ifdef MFEM_USE_MPI
-#
-TEST_CASE("ParGridFunction in Serial",
-          "[ParGridFunction]"
-          "[Parallel]")
+
+TEST_CASE("ParGridFunction in Serial", "[ParGridFunction][Parallel]")
 {
-   int num_procs;
-   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   const int num_procs = Mpi::WorldSize();
+   const int my_rank = Mpi::WorldRank();
 
-   int my_rank;
-   MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
+   const int order = 1;
+   const int save_rank = 0;
 
-   int order = 1;
-   int save_rank = 0;
-   int n = 2 * num_procs;
-
-   FunctionCoefficient squaredFC(squared);
+   const int n = 2 * num_procs;
    Mesh mesh = Mesh::MakeCartesian2D(n, n, Element::QUADRILATERAL);
    mesh.SetCurvature(2);
-   H1_FECollection fec(order, mesh.Dimension());
-
-   double ser_l2_err = 0.0;
-   FiniteElementSpace fespace(&mesh, &fec);
-   GridFunction x(&fespace);
-   x.ProjectCoefficient(squaredFC);
-   ser_l2_err = x.ComputeL2Error(squaredFC);
 
    // Define a parallel mesh by a partitioning of the serial mesh.
    ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   H1_FECollection fec(order, mesh.Dimension());
    ParFiniteElementSpace pfespace(&pmesh, &fec);
    ParGridFunction px(&pfespace);
-   px.ProjectCoefficient(squaredFC);
+   px.Randomize(1);
+   // Ensure that the L-DOFs are set consistently on all ranks
+   px.SetTrueVector();
+   px.SetFromTrueVector();
+
+   ConstantCoefficient zero(0.0);
+   const double l2_norm = px.ComputeL2Error(zero);
 
    // Get the ParMesh and ParGridFunction on 1 of the mpi ranks. Check the
    // L2 error on that rank and save gridfunction.
    Mesh par_to_ser_mesh = pmesh.GetSerialMesh(save_rank);
-   GridFunction x_par_to_ser = px.GetSerialGridFunction(save_rank,
-                                                        par_to_ser_mesh);
-   double par_to_ser_l2_err = 0.0;
+   GridFunction x_par_to_ser = px.GetSerialGridFunction(
+                                  save_rank, par_to_ser_mesh);
    if (my_rank == save_rank)
    {
-      par_to_ser_l2_err = x_par_to_ser.ComputeL2Error(squaredFC);
-      REQUIRE(std::fabs(par_to_ser_l2_err-ser_l2_err) == MFEM_Approx(0.0));
-
-      // Save
+      const double par_to_ser_l2_norm = x_par_to_ser.ComputeL2Error(zero);
+      REQUIRE(par_to_ser_l2_norm == MFEM_Approx(l2_norm));
+      // Save to disk
       par_to_ser_mesh.Save("parallel_in_serial.mesh");
    }
-   MPI_Barrier(MPI_COMM_WORLD);
+
+   {
+      FiniteElementSpace &fes = *x_par_to_ser.FESpace();
+      GridFunction x_par_to_ser_2 = px.GetSerialGridFunction(save_rank, fes);
+      x_par_to_ser_2 -= x_par_to_ser;
+      REQUIRE(x_par_to_ser_2.Normlinf() == MFEM_Approx(0.0));
+   }
 
    // Save the mesh and then load the saved mesh and gridfunction, and check
    // the L2 error on all ranks.
    px.SaveAsSerial("parallel_in_serial.gf", 16, save_rank);
-   Mesh par_to_ser_mesh_read = Mesh("parallel_in_serial.mesh");
-   named_ifgzstream gfstream("parallel_in_serial.gf");
-   GridFunction x_par_to_ser_read = GridFunction(&par_to_ser_mesh_read,
-                                                 gfstream);
-   double par_to_ser_l2_read_err = x_par_to_ser_read.ComputeL2Error(squaredFC);
-   REQUIRE(std::fabs(par_to_ser_l2_read_err-ser_l2_err) == MFEM_Approx(0.0));
-}
-#endif // MFEM_USE_MPI
 
-} // namespace pgridfunc_save_in_serial
+   if (my_rank == save_rank)
+   {
+      Mesh par_to_ser_mesh_read = Mesh("parallel_in_serial.mesh");
+      named_ifgzstream gfstream("parallel_in_serial.gf");
+      GridFunction x_par_to_ser_read(&par_to_ser_mesh_read, gfstream);
+      const double par_to_ser_l2_read_norm = x_par_to_ser_read.ComputeL2Error(zero);
+
+      REQUIRE(par_to_ser_l2_read_norm == MFEM_Approx(l2_norm));
+   }
+
+   if (my_rank == save_rank)
+   {
+      // Clean up
+      REQUIRE(std::remove("parallel_in_serial.mesh") == 0);
+      REQUIRE(std::remove("parallel_in_serial.gf") == 0);
+   }
+}
+
+#endif // MFEM_USE_MPI

tests/unit/general/test_array.cpp

@@ -14,15 +14,23 @@
 
 using namespace mfem;
 
-TEST_CASE("Array init-list construction", "[Array]")
+TEST_CASE("Array init-list and C-style array constructors", "[Array]")
 {
    int ContigData[6] = {6, 5, 4, 3, 2, 1};
+   // Pointer and size construcor
    Array<int> a(ContigData, 6);
-   Array<int> b{6, 5, 4, 3, 2, 1};
+   // Braced-list constructor
+   Array<int> b({6, 5, 4, 3, 2, 1});
+   // Statically sized C-style array constructor
+   Array<int> c(ContigData);
+   // Convertible type constructor
+   Array<int> d({6.0, 5.0, 4.0, 3.0, 2.0, 1.0});
 
    for (int i = 0; i < a.Size(); i++)
    {
       REQUIRE(a[i] == b[i]);
+      REQUIRE(a[i] == c[i]);
+      REQUIRE(a[i] == d[i]);
    }
 }
 

tests/unit/linalg/test_operator.cpp

@@ -104,3 +104,30 @@ TEST_CASE("ConstrainedOperator", "[ConstrainedOperator][Operator]")
    REQUIRE(constrained_mult_application(A, list, x, y_true_zero_transpose, true,
                                         Operator::DiagonalPolicy::DIAG_ZERO) == MFEM_Approx(0.0));
 }
+
+TEST_CASE("Sum and product operators", "[Operator]")
+{
+   const int n = 1;
+   IdentityOperator op_1(n);
+   auto op_2 = MakeOwning<IdentityOperator>(n);
+
+   // op_1 will not be owned, op_2 will be owned
+   SumOperator sum(NonOwning(&op_1), 1.0, op_2, 1.0);
+   ProductOperator product(NonOwning(&op_1), op_2);
+   // Note: it is not a problem for triple to own op_2 'twice'
+   TripleProductOperator triple(NonOwning(&op_1), op_2, op_2);
+
+   // Even though op_2 is reset here, it remains valid in each of the operators
+   op_2.Reset();
+
+   Vector x(n), y(n);
+   x = 1.0;
+   sum.Mult(x, y);
+   REQUIRE(y[0] == 2.0);
+
+   product.Mult(x, y);
+   REQUIRE(y[0] == 1.0);
+
+   triple.Mult(x, y);
+   REQUIRE(y[0] == 1.0);
+}

tests/unit/linalg/test_vector.cpp

@@ -15,15 +15,23 @@
 
 using namespace mfem;
 
-TEST_CASE("Vector init-list construction", "[Vector]")
+TEST_CASE("Vector init-list and C-style array constructors", "[Vector]")
 {
    real_t ContigData[6] = {6.0, 5.0, 4.0, 3.0, 2.0, 1.0};
+   // Point and size constructor
    Vector a(ContigData, 6);
+   // Braced-list constructor
    Vector b({6.0, 5.0, 4.0, 3.0, 2.0, 1.0});
+   // Statically sized C-style array constructor
+   Vector c(ContigData);
+   // Convertible type constructor
+   Vector d({6, 5, 4, 3, 2, 1});
 
    for (int i = 0; i < a.Size(); i++)
    {
-      REQUIRE(a(i) == b(i));
+      REQUIRE(a[i] == b[i]);
+      REQUIRE(a[i] == c[i]);
+      REQUIRE(a[i] == d[i]);
    }
 }
 

tests/unit/mesh/test_mesh.cpp

@@ -180,3 +180,88 @@ TEST_CASE("MakeSimplicial", "[Mesh]")
    // on the original mesh, but it doesn't happen for these test cases.
    REQUIRE(simplex_mesh.GetNE() == orig_mesh.GetNE()*factor);
 }
+
+TEST_CASE("MakeNurbs", "[Mesh]")
+{
+   Array<real_t> intervals_array({1, 1, 1});
+   Vector intervals(intervals_array.GetData(), intervals_array.Size());
+   Array<int> continuity({-1, 1, 1, -1});
+   {
+      const KnotVector kv(2, intervals, continuity);
+      REQUIRE(kv.GetNE() == 3);
+      REQUIRE(kv.GetNCP() == 5);
+      REQUIRE(kv.GetOrder() == 2);
+      REQUIRE(kv.Size() == 8);
+   }
+   {
+      const KnotVector kv(3, intervals, continuity);
+      REQUIRE(kv.GetNE() == 3);
+      REQUIRE(kv.GetNCP() == 8);
+      REQUIRE(kv.GetOrder() == 3);
+      REQUIRE(kv.Size() == 12);
+   }
+   const KnotVector kv(2, intervals, continuity);
+   Array<real_t> grev_pts({0.0, 1.0/6.0, 0.5, 5.0/6.0, 1.0});
+
+   Array<NURBSPatch *> patches;
+
+   // Will build and test on multiple NURBS meshes. Cleans up and
+   // resets the patches array, which is assumed to be initially
+   // populated for the particular test case.
+   const auto test_nurbs_extension = [&](Mesh& patch_topology)
+   {
+      NURBSExtension ne(&patch_topology, patches);
+      Mesh mesh(ne);
+      GridFunction *nodes = mesh.GetNodes();
+      REQUIRE(nodes != NULL);
+      FiniteElementSpace *fe = nodes->FESpace();
+      REQUIRE(fe != NULL);
+      SparseMatrix p(fe->GetNDofs(), fe->GetNDofs());
+      SparseMatrix r(fe->GetNDofs(), fe->GetNDofs());
+      for (int i = 0; i < fe->GetNDofs(); ++i)
+      {
+         p.Add(i, i, 1);
+         r.Add(i, i, 1);
+      }
+      p.Finalize();
+      r.Finalize();
+      fe->SetProlongation(p);
+      fe->SetRestriction(r);
+      for (int i=0; i<patches.Size(); i++) { delete patches[i]; }
+      patches.SetSize(0);
+   };
+
+   // Bi-variate 2D test:
+   Array<real_t> pts_2d(3 * kv.GetNCP() * kv.GetNCP());
+   int count = 0;
+   for (int j = 0; j < kv.GetNCP(); ++j)
+      for (int i = 0; i < kv.GetNCP(); ++i)
+      {
+         pts_2d[count + 0] = grev_pts[i];
+         pts_2d[count + 1] = grev_pts[j];
+         pts_2d[count + 2] = 1;
+         count += 3;
+      }
+   patches.Append(new NURBSPatch(&kv, &kv, 3, pts_2d.GetData()));
+   Mesh patch_topology_2d =
+      Mesh::MakeCartesian2D(1, 1, Element::Type::QUADRILATERAL);
+   test_nurbs_extension(patch_topology_2d);
+
+   // Tri-variate 3D test:
+   Array<real_t> pts_3d(4 * kv.GetNCP() * kv.GetNCP() * kv.GetNCP());
+   count = 0;
+   for (int k = 0; k < kv.GetNCP(); ++k)
+      for (int j = 0; j < kv.GetNCP(); ++j)
+         for (int i = 0; i < kv.GetNCP(); ++i)
+         {
+            pts_3d[count + 0] = grev_pts[i];
+            pts_3d[count + 1] = grev_pts[j];
+            pts_3d[count + 2] = grev_pts[k];
+            pts_3d[count + 3] = 1;
+            count += 4;
+         }
+   patches.Append(new NURBSPatch(&kv, &kv, &kv, 4, pts_3d.GetData()));
+   Mesh patch_topology_3d =
+      Mesh::MakeCartesian3D(1, 1, 1, Element::Type::HEXAHEDRON);
+   test_nurbs_extension(patch_topology_3d);
+}