make style

fem/bilinearform.cpp

@@ -100,6 +100,7 @@ BilinearForm::BilinearForm (FiniteElementSpace * f, BilinearForm * bf, int ps)
 
    // Copy the pointers to the integrators
    domain_integs = bf->domain_integs;
+   domain_integs_marker = bf->domain_integs_marker;
 
    boundary_integs = bf->boundary_integs;
    boundary_integs_marker = bf->boundary_integs_marker;
@@ -112,6 +113,138 @@ BilinearForm::BilinearForm (FiniteElementSpace * f, BilinearForm * bf, int ps)
    AllocMat();
 }
 
+BilinearForm::BilinearForm(BilinearForm &&other)
+   : Matrix(other.fes->GetVSize()), mat(other.mat), mat_e(other.mat_e),
+     fes(other.fes), assembly(other.assembly), batch(other.batch),
+     ext(other.ext), sequence(other.sequence), extern_bfs(other.extern_bfs),
+     element_matrices(other.element_matrices), static_cond(other.static_cond),
+     hybridization(other.hybridization), diag_policy(other.diag_policy),
+     precompute_sparsity(other.precompute_sparsity)
+{
+   // We swap stored integrators and markers with the moved nonlinear form
+   mfem::Swap(domain_integs, other.domain_integs);
+   mfem::Swap(domain_integs_marker, other.domain_integs_marker);
+   mfem::Swap(boundary_integs, other.boundary_integs);
+   mfem::Swap(boundary_integs_marker, other.boundary_integs_marker);
+   mfem::Swap(interior_face_integs, other.interior_face_integs);
+   mfem::Swap(boundary_face_integs, other.boundary_face_integs);
+   mfem::Swap(boundary_face_integs_marker, other.boundary_face_integs_marker);
+
+   /// Leave the moved nonlinear form in a state as if it was just constructed
+   /// with fes
+   other.sequence = fes->GetSequence();
+   other.mat = nullptr;
+   other.mat_e = nullptr;
+   other.extern_bfs = 0;
+   other.element_matrices = nullptr;
+   other.static_cond = nullptr;
+   other.hybridization = nullptr;
+   other.precompute_sparsity = 0;
+   other.diag_policy = DIAG_KEEP;
+
+   other.assembly = AssemblyLevel::LEGACY;
+   other.batch = 1;
+   other.ext = nullptr;
+}
+
+BilinearForm& BilinearForm::operator=(BilinearForm &&other)
+{
+   if (this != &other)
+   {
+      /// Cleanup current bilinear form first
+      delete mat_e;
+      delete mat;
+      delete element_matrices;
+      delete static_cond;
+      delete hybridization;
+      for (int i = 0; i < domain_integs.Size(); i++) { delete domain_integs[i]; }
+      for (int i = 0; i < boundary_integs.Size(); i++) { delete boundary_integs[i]; }
+      for (int i = 0; i < interior_face_integs.Size(); i++) { delete interior_face_integs[i]; }
+      for (int i = 0; i < boundary_face_integs.Size(); i++) { delete boundary_face_integs[i]; }
+      delete ext;
+
+      /// Null out all our integs and set size of their arrays to zero
+      for (int k = 0; k < domain_integs.Size(); k++)
+      {
+         domain_integs[k] = nullptr;
+      }
+      domain_integs.SetSize(0);
+      for (int k = 0; k < boundary_integs.Size(); k++)
+      {
+         boundary_integs[k] = nullptr;
+      }
+      boundary_integs.SetSize(0);
+      for (int k = 0; k < interior_face_integs.Size(); k++)
+      {
+         interior_face_integs[k] = nullptr;
+      }
+      interior_face_integs.SetSize(0);
+      for (int k = 0; k < boundary_face_integs.Size(); ++k)
+      {
+         boundary_face_integs[k] = nullptr;
+      }
+      boundary_face_integs.SetSize(0);
+
+      /// Null out all our markers and set size of their arrays to zero
+      for (int k = 0; k < domain_integs_marker.Size(); ++k)
+      {
+         domain_integs_marker[k] = nullptr;
+      }
+      domain_integs_marker.SetSize(0);
+      for (int k = 0; k < boundary_integs_marker.Size(); ++k)
+      {
+         boundary_integs_marker[k] = nullptr;
+      }
+      boundary_integs_marker.SetSize(0);
+      for (int k = 0; k < boundary_face_integs_marker.Size(); ++k)
+      {
+         boundary_face_integs_marker[k] = nullptr;
+      }
+      boundary_face_integs_marker.SetSize(0);
+
+      /// Now steal data from other bilinear form leaving it in a state as if
+      /// it was just constructed with fes
+      Matrix::operator=(std::move(other));
+
+      mat = other.mat;
+      other.mat = nullptr;
+      mat_e = other.mat_e;
+      other.mat_e = nullptr;
+      fes = other.fes;
+      assembly = other.assembly;
+      other.assembly = AssemblyLevel::LEGACY;
+      batch = other.batch;
+      other.batch = 1;
+      ext = other.ext;
+      other.ext = nullptr;
+      sequence = other.sequence;
+      other.sequence = fes->GetSequence();
+      extern_bfs = other.extern_bfs;
+      other.extern_bfs = 0;
+
+      // Swap our empty integ and marker arrays with the moved bilinear form
+      mfem::Swap(domain_integs, other.domain_integs);
+      mfem::Swap(domain_integs_marker, other.domain_integs_marker);
+      mfem::Swap(boundary_integs, other.boundary_integs);
+      mfem::Swap(boundary_integs_marker, other.boundary_integs_marker);
+      mfem::Swap(interior_face_integs, other.interior_face_integs);
+      mfem::Swap(boundary_face_integs, other.boundary_face_integs);
+      mfem::Swap(boundary_face_integs_marker, other.boundary_face_integs_marker);
+
+      element_matrices = other.element_matrices;
+      other.element_matrices = nullptr;
+      static_cond = other.static_cond;
+      other.static_cond = nullptr;
+      hybridization = other.hybridization;
+      other.hybridization = nullptr;
+      diag_policy = other.diag_policy;
+      other.diag_policy = DIAG_KEEP;
+      precompute_sparsity = other.precompute_sparsity;
+      other.precompute_sparsity = 0;
+   }
+   return *this;
+}
+
 void BilinearForm::SetAssemblyLevel(AssemblyLevel assembly_level)
 {
    if (ext)

fem/bilinearform.hpp

@@ -90,23 +90,23 @@ protected:
    int extern_bfs;
 
    /// Set of Domain Integrators to be applied.
-   Array<BilinearFormIntegrator*> domain_integs;
+   Array<BilinearFormIntegrator*> domain_integs; // owned
    /// Element attribute marker (should be of length mesh->attributes.Max() or
    /// 0 if mesh->attributes is empty)
    /// Includes all by default.
    /// 0 - ignore attribute
    /// 1 - include attribute
-   Array<Array<int>*>             domain_integs_marker;
+   Array<Array<int>*>             domain_integs_marker; // not owned
 
    /// Set of Boundary Integrators to be applied.
-   Array<BilinearFormIntegrator*> boundary_integs;
+   Array<BilinearFormIntegrator*> boundary_integs; // owned
    Array<Array<int>*> boundary_integs_marker; ///< Entries are not owned.
 
    /// Set of interior face Integrators to be applied.
-   Array<BilinearFormIntegrator*> interior_face_integs;
+   Array<BilinearFormIntegrator*> interior_face_integs; // owned
 
    /// Set of boundary face Integrators to be applied.
-   Array<BilinearFormIntegrator*> boundary_face_integs;
+   Array<BilinearFormIntegrator*> boundary_face_integs; // owned
    Array<Array<int>*> boundary_face_integs_marker; ///< Entries are not owned.
 
    DenseMatrix elemmat;
@@ -141,13 +141,6 @@ protected:
       ext = NULL;
    }
 
-private:
-   /// Copy construction is not supported; body is undefined.
-   BilinearForm(const BilinearForm &);
-
-   /// Copy assignment is not supported; body is undefined.
-   BilinearForm &operator=(const BilinearForm &);
-
 public:
    /// Creates bilinear form associated with FE space @a *f.
    /** The pointer @a f is not owned by the newly constructed object. */
@@ -165,6 +158,20 @@ public:
        #precompute_sparsity, see UsePrecomputedSparsity() for details. */
    BilinearForm(FiniteElementSpace *f, BilinearForm *bf, int ps = 0);
 
+   /// Explicitly prohibit copy construction/assignment of BilinearForms
+   BilinearForm(const BilinearForm &) = delete;
+   BilinearForm &operator=(const BilinearForm &) = delete;
+
+   /// Move constructor for BilinearForm.
+   /** This constructor "steals" the owned data members from the @a other
+       BilinearForm. */
+   BilinearForm(BilinearForm &&other);
+
+   /// Move assignment operator for BilinearForm
+   /** This assignment first frees all owned data, then "steals" the owned data
+       members from the @a other BilinearForm. */
+   BilinearForm& operator=(BilinearForm &&other);
+
    /// Get the size of the BilinearForm as a square matrix.
    int Size() const { return height; }
 

fem/linearform.cpp

@@ -26,15 +26,134 @@ LinearForm::LinearForm(FiniteElementSpace *f, LinearForm *lf)
    ext = nullptr;
    extern_lfs = 1;
 
-   // Copy the pointers to the integrators
+   // Copy the pointers to the integrators and markers
    domain_integs = lf->domain_integs;
+   domain_integs_marker = lf->domain_integs_marker;
 
    domain_delta_integs = lf->domain_delta_integs;
 
    boundary_integs = lf->boundary_integs;
+   boundary_integs_marker = lf->boundary_integs_marker;
 
    boundary_face_integs = lf->boundary_face_integs;
    boundary_face_integs_marker = lf->boundary_face_integs_marker;
+
+   interior_face_integs = lf->interior_face_integs;
+}
+
+LinearForm::LinearForm(LinearForm &&other)
+   : Vector(std::move(other)), fes(other.fes), ext(other.ext),
+     extern_lfs(other.extern_lfs)
+{
+   // Linear forms are stored on the device
+   UseDevice(true);
+
+   // We swap stored integrators and markers with the moved LinearForm
+   mfem::Swap(domain_integs, other.domain_integs);
+   mfem::Swap(domain_integs_marker, other.domain_integs_marker);
+   mfem::Swap(domain_delta_integs, other.domain_delta_integs);
+   mfem::Swap(boundary_integs, other.boundary_integs);
+   mfem::Swap(boundary_integs_marker, other.boundary_integs_marker);
+   mfem::Swap(boundary_face_integs, other.boundary_face_integs);
+   mfem::Swap(boundary_face_integs_marker, other.boundary_face_integs_marker);
+   mfem::Swap(interior_face_integs, other.interior_face_integs);
+
+   other.fes = nullptr;
+
+   other.ext = nullptr;
+
+   // moved LinearForm now set to owns its integrators, though none are stored
+   // since we swapped our newly constructed empty Arrays with its Arrays
+   other.extern_lfs = 0;
+}
+
+LinearForm& LinearForm::operator=(LinearForm &&other)
+{
+   if (this != &other)
+   {
+      Vector::operator=(std::move(other));
+
+      /// If we own our integrators, delete them
+      if (extern_lfs == 0)
+      {
+         for (int k = 0; k < domain_delta_integs.Size(); k++)
+         { delete domain_delta_integs[k]; }
+         for (int k = 0; k < domain_integs.Size(); k++) { delete domain_integs[k]; }
+         for (int k = 0; k < boundary_integs.Size(); k++) { delete boundary_integs[k]; }
+         for (int k = 0; k < boundary_face_integs.Size(); k++)
+         { delete boundary_face_integs[k]; }
+         for (int k = 0; k < interior_face_integs.Size(); k++)
+         { delete interior_face_integs[k]; }
+      }
+
+      /// Null out all integs
+      for (int k = 0; k < domain_delta_integs.Size(); k++)
+      {
+         domain_delta_integs[k] = nullptr;
+      }
+      domain_delta_integs.SetSize(0);
+      for (int k = 0; k < domain_integs.Size(); k++)
+      {
+         domain_integs[k] = nullptr;
+      }
+      domain_integs.SetSize(0);
+      for (int k = 0; k < boundary_integs.Size(); k++)
+      {
+         boundary_integs[k] = nullptr;
+      }
+      boundary_integs.SetSize(0);
+      for (int k = 0; k < boundary_face_integs.Size(); k++)
+      {
+         boundary_face_integs[k] = nullptr;
+      }
+      boundary_face_integs.SetSize(0);
+      for (int k = 0; k < interior_face_integs.Size(); k++)
+      {
+         interior_face_integs[k] = nullptr;
+      }
+      interior_face_integs.SetSize(0);
+
+      /// Null out all markers
+      for (int k = 0; k < domain_integs_marker.Size(); ++k)
+      {
+         domain_integs_marker[k] = nullptr;
+      }
+      domain_integs_marker.SetSize(0);
+      for (int k = 0; k < boundary_integs_marker.Size(); ++k)
+      {
+         boundary_integs_marker[k] = nullptr;
+      }
+      boundary_integs_marker.SetSize(0);
+      for (int k = 0; k < boundary_face_integs_marker.Size(); ++k)
+      {
+         boundary_face_integs_marker[k] = nullptr;
+      }
+      boundary_face_integs_marker.SetSize(0);
+
+      Vector::operator=(std::move(other));
+
+      // swap stored integrators and markers with the moved LinearForm
+      mfem::Swap(domain_integs, other.domain_integs);
+      mfem::Swap(domain_integs_marker, other.domain_integs_marker);
+      mfem::Swap(domain_delta_integs, other.domain_delta_integs);
+      mfem::Swap(boundary_integs, other.boundary_integs);
+      mfem::Swap(boundary_integs_marker, other.boundary_integs_marker);
+      mfem::Swap(boundary_face_integs, other.boundary_face_integs);
+      mfem::Swap(boundary_face_integs_marker, other.boundary_face_integs_marker);
+      mfem::Swap(interior_face_integs, other.interior_face_integs);
+
+      fes = other.fes;
+      other.fes = nullptr;
+
+      ext = other.ext;
+      other.ext = nullptr;
+
+      extern_lfs = other.extern_lfs;
+      // moved LinearForm now set to owns its integrators, though none are stored
+      // since we swapped our empty Arrays with its Arrays
+      other.extern_lfs = 0;
+   }
+   return *this;
 }
 
 void LinearForm::AddDomainIntegrator(LinearFormIntegrator *lfi)

fem/linearform.hpp

@@ -74,10 +74,6 @@ protected:
    /// Force (re)computation of delta locations.
    void ResetDeltaLocations() { domain_delta_integs_elem_id.SetSize(0); }
 
-private:
-   /// Copy construction is not supported; body is undefined.
-   LinearForm(const LinearForm &);
-
 public:
    /// Creates linear form associated with FE space @a *f.
    /** The pointer @a f is not owned by the newly constructed object. */
@@ -107,6 +103,9 @@ public:
    LinearForm(FiniteElementSpace *f, double *data) : Vector(data, f->GetVSize())
    { fes = f; ext = nullptr; extern_lfs = 0; }
 
+   /// Explicitly prohibit copy construction of LinearForm
+   LinearForm(const LinearForm &other) = delete;
+
    /// Copy assignment. Only the data of the base class Vector is copied.
    /** It is assumed that this object and @a rhs use FiniteElementSpace%s that
        have the same size.
@@ -116,6 +115,16 @@ public:
    LinearForm &operator=(const LinearForm &rhs)
    { return operator=((const Vector &)rhs); }
 
+   /// Move constructor for LinearForm.
+   /** This constructor "steals" the owned data members from the @a other
+       LinearForm. */
+   LinearForm(LinearForm &&other);
+
+   /// Move assignment operator for LinearForm
+   /** This assignment first frees all owned data, then "steals" the owned data
+       members from the @a other LinearForm. */
+   LinearForm& operator=(LinearForm &&other);
+
    /// (DEPRECATED) Return the FE space associated with the LinearForm.
    /** @deprecated Use FESpace() instead. */
    MFEM_DEPRECATED FiniteElementSpace *GetFES() { return fes; }

fem/nonlinearform.cpp

@@ -15,6 +15,92 @@
 namespace mfem
 {
 
+NonlinearForm::NonlinearForm(NonlinearForm &&other)
+   : Operator(other.fes->GetTrueVSize()), assembly(other.assembly),
+     ext(other.ext), fes(other.fes), Grad(other.Grad), cGrad(other.cGrad),
+     sequence(other.fes->GetSequence()), P(other.fes->GetProlongationMatrix()),
+     cP(dynamic_cast<const SparseMatrix*>(P))
+{
+   // We swap stored integrators and markers with the moved nonlinear form
+   mfem::Swap(domain_integs, other.domain_integs);
+   mfem::Swap(domain_integs_marker, other.domain_integs_marker);
+   mfem::Swap(interior_face_integs, other.interior_face_integs);
+   mfem::Swap(boundary_face_integs, other.boundary_face_integs);
+   mfem::Swap(boundary_face_integs_marker, other.boundary_face_integs_marker);
+
+   /// Leave the moved nonlinear form in a state as if it was just constructed
+   /// with fes
+   other.ext = nullptr;
+   other.cGrad = nullptr;
+   other.Grad = nullptr;
+   other.assembly = AssemblyLevel::LEGACY;
+}
+
+NonlinearForm& NonlinearForm::operator=(NonlinearForm &&other)
+{
+   if (this != &other)
+   {
+      /// Cleanup current nonlinear form first
+      delete cGrad;
+      delete Grad;
+      for (int i = 0; i < domain_integs.Size(); i++) { delete domain_integs[i]; }
+      for (int i = 0; i < interior_face_integs.Size(); i++) { delete interior_face_integs[i]; }
+      for (int i = 0; i < boundary_face_integs.Size(); i++) { delete boundary_face_integs[i]; }
+      delete ext;
+
+      /// Null out all our integs and set size of their arrays to zero
+      for (int k = 0; k < domain_integs.Size(); k++)
+      {
+         domain_integs[k] = nullptr;
+      }
+      domain_integs.SetSize(0);
+      for (int k = 0; k < boundary_face_integs.Size(); k++)
+      {
+         boundary_face_integs[k] = nullptr;
+      }
+      boundary_face_integs.SetSize(0);
+      for (int k = 0; k < interior_face_integs.Size(); k++)
+      {
+         interior_face_integs[k] = nullptr;
+      }
+      interior_face_integs.SetSize(0);
+
+      /// Null out all our markers and set size of their arrays to zero
+      for (int k = 0; k < domain_integs_marker.Size(); ++k)
+      {
+         domain_integs_marker[k] = nullptr;
+      }
+      domain_integs_marker.SetSize(0);
+      for (int k = 0; k < boundary_face_integs_marker.Size(); ++k)
+      {
+         boundary_face_integs_marker[k] = nullptr;
+      }
+      boundary_face_integs_marker.SetSize(0);
+
+      /// Now steal data from other nonlinear form leaving it in a state as if
+      /// it was just constructed with fes
+      Operator::operator=(std::move(other));
+
+      assembly = other.assembly;
+      other.assembly = AssemblyLevel::LEGACY;
+      Grad = other.Grad;
+      other.Grad = nullptr;
+      cGrad = other.cGrad;
+      other.cGrad = nullptr;
+
+      // Swap our empty integ and marker arrays with the moved nonlinear form
+      mfem::Swap(domain_integs, other.domain_integs);
+      mfem::Swap(domain_integs_marker, other.domain_integs_marker);
+      mfem::Swap(interior_face_integs, other.interior_face_integs);
+      mfem::Swap(boundary_face_integs, other.boundary_face_integs);
+      mfem::Swap(boundary_face_integs_marker, other.boundary_face_integs_marker);
+
+      ext = other.ext;
+      other.ext = nullptr;
+   }
+   return *this;
+}
+
 void NonlinearForm::SetAssemblyLevel(AssemblyLevel assembly_level)
 {
    if (ext)
@@ -38,6 +124,38 @@ void NonlinearForm::SetAssemblyLevel(AssemblyLevel assembly_level)
    }
 }
 
+void NonlinearForm::AddDomainIntegrator(NonlinearFormIntegrator *nlfi)
+{
+   domain_integs.Append(nlfi);
+   domain_integs_marker.Append(nullptr); // null marker means apply everywhere
+}
+
+void NonlinearForm::AddDomainIntegrator(NonlinearFormIntegrator *nlfi,
+                                        Array<int> &elem_marker)
+{
+   domain_integs.Append(nlfi);
+   domain_integs_marker.Append(&elem_marker);
+}
+
+void NonlinearForm::AddInteriorFaceIntegrator(NonlinearFormIntegrator *nlfi)
+{
+   interior_face_integs.Append(nlfi);
+}
+
+void NonlinearForm::AddBdrFaceIntegrator(NonlinearFormIntegrator *nlfi)
+{
+   boundary_face_integs.Append(nlfi);
+   // null marker means apply everywhere
+   boundary_face_integs_marker.Append(nullptr);
+}
+
+void NonlinearForm::AddBdrFaceIntegrator(NonlinearFormIntegrator *nlfi,
+                                         Array<int> &bdr_marker)
+{
+   boundary_face_integs.Append(nlfi);
+   boundary_face_integs_marker.Append(&bdr_marker);
+}
+
 void NonlinearForm::SetEssentialBC(const Array<int> &bdr_attr_is_ess,
                                    Vector *rhs)
 {
@@ -87,8 +205,10 @@ double NonlinearForm::GetGridFunctionEnergy(const Vector &x) const
 {
    if (ext)
    {
-      MFEM_VERIFY(!fnfi.Size(), "Interior faces terms not yet implemented!");
-      MFEM_VERIFY(!bfnfi.Size(), "Boundary face terms not yet implemented!");
+      MFEM_VERIFY(!interior_face_integs.Size(),
+                  "Interior faces terms not yet implemented!");
+      MFEM_VERIFY(!boundary_face_integs.Size(),
+                  "Boundary face terms not yet implemented!");
       return ext->GetGridFunctionEnergy(x);
    }
 
@@ -97,30 +217,47 @@ double NonlinearForm::GetGridFunctionEnergy(const Vector &x) const
    const FiniteElement *fe;
    ElementTransformation *T;
    DofTransformation *doftrans;
+   Mesh *mesh = fes->GetMesh();
    double energy = 0.0;
 
-   if (dnfi.Size())
+   if (domain_integs.Size())
    {
+      for (int k = 0; k < domain_integs.Size(); k++)
+      {
+         if (domain_integs_marker[k] != nullptr)
+         {
+            MFEM_VERIFY(mesh->attributes.Size() ==
+                        domain_integs_marker[k]->Size(),
+                        "invalid element marker for domain integrator #"
+                        << k << ", counting from zero");
+         }
+      }
+
       for (int i = 0; i < fes->GetNE(); i++)
       {
+         int elem_attr = mesh->GetAttribute(i);
          fe = fes->GetFE(i);
          doftrans = fes->GetElementVDofs(i, vdofs);
          T = fes->GetElementTransformation(i);
          x.GetSubVector(vdofs, el_x);
          if (doftrans) {doftrans->InvTransformPrimal(el_x); }
-         for (int k = 0; k < dnfi.Size(); k++)
+         for (int k = 0; k < domain_integs.Size(); k++)
          {
-            energy += dnfi[k]->GetElementEnergy(*fe, *T, el_x);
+            if (domain_integs_marker[k] == nullptr ||
+                (*(domain_integs_marker[k]))[elem_attr-1] == 1)
+            {
+               energy += domain_integs[k]->GetElementEnergy(*fe, *T, el_x);
+            }
          }
       }
    }
 
-   if (fnfi.Size())
+   if (interior_face_integs.Size())
    {
       MFEM_ABORT("TODO: add energy contribution from interior face terms");
    }
 
-   if (bfnfi.Size())
+   if (boundary_face_integs.Size())
    {
       MFEM_ABORT("TODO: add energy contribution from boundary face terms");
    }
@@ -173,25 +310,41 @@ void NonlinearForm::Mult(const Vector &x, Vector &y) const
 
    py = 0.0;
 
-   if (dnfi.Size())
+   if (domain_integs.Size())
    {
+      for (int k = 0; k < domain_integs.Size(); k++)
+      {
+         if (domain_integs_marker[k] != nullptr)
+         {
+            MFEM_VERIFY(mesh->attributes.Size() ==
+                        domain_integs_marker[k]->Size(),
+                        "invalid element marker for domain integrator #"
+                        << k << ", counting from zero");
+         }
+      }
+
       for (int i = 0; i < fes->GetNE(); i++)
       {
+         int elem_attr = mesh->GetAttribute(i);
          fe = fes->GetFE(i);
          doftrans = fes->GetElementVDofs(i, vdofs);
          T = fes->GetElementTransformation(i);
          px.GetSubVector(vdofs, el_x);
          if (doftrans) {doftrans->InvTransformPrimal(el_x); }
-         for (int k = 0; k < dnfi.Size(); k++)
+         for (int k = 0; k < domain_integs.Size(); k++)
          {
-            dnfi[k]->AssembleElementVector(*fe, *T, el_x, el_y);
-            if (doftrans) {doftrans->TransformDual(el_y); }
-            py.AddElementVector(vdofs, el_y);
+            if (domain_integs_marker[k] == nullptr ||
+                (*(domain_integs_marker[k]))[elem_attr-1] == 1)
+            {
+               domain_integs[k]->AssembleElementVector(*fe, *T, el_x, el_y);
+               if (doftrans) {doftrans->TransformDual(el_y); }
+               py.AddElementVector(vdofs, el_y);
+            }
          }
       }
    }
 
-   if (fnfi.Size())
+   if (interior_face_integs.Size())
    {
       FaceElementTransformations *tr;
       const FiniteElement *fe1, *fe2;
@@ -211,16 +364,16 @@ void NonlinearForm::Mult(const Vector &x, Vector &y) const
             fe1 = fes->GetFE(tr->Elem1No);
             fe2 = fes->GetFE(tr->Elem2No);
 
-            for (int k = 0; k < fnfi.Size(); k++)
+            for (int k = 0; k < interior_face_integs.Size(); k++)
             {
-               fnfi[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
+               interior_face_integs[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
                py.AddElementVector(vdofs, el_y);
             }
          }
       }
    }
 
-   if (bfnfi.Size())
+   if (boundary_face_integs.Size())
    {
       FaceElementTransformations *tr;
       const FiniteElement *fe1, *fe2;
@@ -229,14 +382,14 @@ void NonlinearForm::Mult(const Vector &x, Vector &y) const
       Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
                                  mesh->bdr_attributes.Max() : 0);
       bdr_attr_marker = 0;
-      for (int k = 0; k < bfnfi.Size(); k++)
+      for (int k = 0; k < boundary_face_integs.Size(); k++)
       {
-         if (bfnfi_marker[k] == NULL)
+         if (boundary_face_integs_marker[k] == NULL)
          {
             bdr_attr_marker = 1;
             break;
          }
-         Array<int> &bdr_marker = *bfnfi_marker[k];
+         Array<int> &bdr_marker = *boundary_face_integs_marker[k];
          MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
                      "invalid boundary marker for boundary face integrator #"
                      << k << ", counting from zero");
@@ -262,12 +415,12 @@ void NonlinearForm::Mult(const Vector &x, Vector &y) const
             // but we can't dereference a NULL pointer, and we don't want to
             // actually make a fake element.
             fe2 = fe1;
-            for (int k = 0; k < bfnfi.Size(); k++)
+            for (int k = 0; k < boundary_face_integs.Size(); k++)
             {
-               if (bfnfi_marker[k] &&
-                   (*bfnfi_marker[k])[bdr_attr-1] == 0) { continue; }
+               if (boundary_face_integs_marker[k] &&
+                   (*boundary_face_integs_marker[k])[bdr_attr-1] == 0) { continue; }
 
-               bfnfi[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
+               boundary_face_integs[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
                py.AddElementVector(vdofs, el_y);
             }
          }
@@ -320,26 +473,42 @@ Operator &NonlinearForm::GetGradient(const Vector &x) const
       *Grad = 0.0;
    }
 
-   if (dnfi.Size())
+   if (domain_integs.Size())
    {
+      for (int k = 0; k < domain_integs.Size(); k++)
+      {
+         if (domain_integs_marker[k] != nullptr)
+         {
+            MFEM_VERIFY(mesh->attributes.Size() ==
+                        domain_integs_marker[k]->Size(),
+                        "invalid element marker for domain integrator #"
+                        << k << ", counting from zero");
+         }
+      }
+
       for (int i = 0; i < fes->GetNE(); i++)
       {
+         int elem_attr = fes->GetMesh()->GetAttribute(i);
          fe = fes->GetFE(i);
          doftrans = fes->GetElementVDofs(i, vdofs);
          T = fes->GetElementTransformation(i);
          px.GetSubVector(vdofs, el_x);
          if (doftrans) {doftrans->InvTransformPrimal(el_x); }
-         for (int k = 0; k < dnfi.Size(); k++)
+         for (int k = 0; k < domain_integs.Size(); k++)
          {
-            dnfi[k]->AssembleElementGrad(*fe, *T, el_x, elmat);
-            if (doftrans) { doftrans->TransformDual(elmat); }
-            Grad->AddSubMatrix(vdofs, vdofs, elmat, skip_zeros);
-            // Grad->AddSubMatrix(vdofs, vdofs, elmat, 1);
+            if (domain_integs_marker[k] == nullptr ||
+                (*(domain_integs_marker[k]))[elem_attr-1] == 1)
+            {
+               domain_integs[k]->AssembleElementGrad(*fe, *T, el_x, elmat);
+               if (doftrans) { doftrans->TransformDual(elmat); }
+               Grad->AddSubMatrix(vdofs, vdofs, elmat, skip_zeros);
+               // Grad->AddSubMatrix(vdofs, vdofs, elmat, 1);
+            }
          }
       }
    }
 
-   if (fnfi.Size())
+   if (interior_face_integs.Size())
    {
       FaceElementTransformations *tr;
       const FiniteElement *fe1, *fe2;
@@ -359,16 +528,16 @@ Operator &NonlinearForm::GetGradient(const Vector &x) const
             fe1 = fes->GetFE(tr->Elem1No);
             fe2 = fes->GetFE(tr->Elem2No);
 
-            for (int k = 0; k < fnfi.Size(); k++)
+            for (int k = 0; k < interior_face_integs.Size(); k++)
             {
-               fnfi[k]->AssembleFaceGrad(*fe1, *fe2, *tr, el_x, elmat);
+               interior_face_integs[k]->AssembleFaceGrad(*fe1, *fe2, *tr, el_x, elmat);
                Grad->AddSubMatrix(vdofs, vdofs, elmat, skip_zeros);
             }
          }
       }
    }
 
-   if (bfnfi.Size())
+   if (boundary_face_integs.Size())
    {
       FaceElementTransformations *tr;
       const FiniteElement *fe1, *fe2;
@@ -377,14 +546,14 @@ Operator &NonlinearForm::GetGradient(const Vector &x) const
       Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
                                  mesh->bdr_attributes.Max() : 0);
       bdr_attr_marker = 0;
-      for (int k = 0; k < bfnfi.Size(); k++)
+      for (int k = 0; k < boundary_face_integs.Size(); k++)
       {
-         if (bfnfi_marker[k] == NULL)
+         if (boundary_face_integs_marker[k] == NULL)
          {
             bdr_attr_marker = 1;
             break;
          }
-         Array<int> &bdr_marker = *bfnfi_marker[k];
+         Array<int> &bdr_marker = *boundary_face_integs_marker[k];
          MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
                      "invalid boundary marker for boundary face integrator #"
                      << k << ", counting from zero");
@@ -410,12 +579,12 @@ Operator &NonlinearForm::GetGradient(const Vector &x) const
             // but we can't dereference a NULL pointer, and we don't want to
             // actually make a fake element.
             fe2 = fe1;
-            for (int k = 0; k < bfnfi.Size(); k++)
+            for (int k = 0; k < boundary_face_integs.Size(); k++)
             {
-               if (bfnfi_marker[k] &&
-                   (*bfnfi_marker[k])[bdr_attr-1] == 0) { continue; }
+               if (boundary_face_integs_marker[k] &&
+                   (*boundary_face_integs_marker[k])[bdr_attr-1] == 0) { continue; }
 
-               bfnfi[k]->AssembleFaceGrad(*fe1, *fe2, *tr, el_x, elmat);
+               boundary_face_integs[k]->AssembleFaceGrad(*fe1, *fe2, *tr, el_x, elmat);
                Grad->AddSubMatrix(vdofs, vdofs, elmat, skip_zeros);
             }
          }
@@ -471,9 +640,9 @@ NonlinearForm::~NonlinearForm()
 {
    delete cGrad;
    delete Grad;
-   for (int i = 0; i <  dnfi.Size(); i++) { delete  dnfi[i]; }
-   for (int i = 0; i <  fnfi.Size(); i++) { delete  fnfi[i]; }
-   for (int i = 0; i < bfnfi.Size(); i++) { delete bfnfi[i]; }
+   for (int i = 0; i < domain_integs.Size(); i++) { delete domain_integs[i]; }
+   for (int i = 0; i < interior_face_integs.Size(); i++) { delete interior_face_integs[i]; }
+   for (int i = 0; i < boundary_face_integs.Size(); i++) { delete boundary_face_integs[i]; }
    delete ext;
 }
 
@@ -564,8 +733,8 @@ BlockNonlinearForm::BlockNonlinearForm(Array<FiniteElementSpace *> &f) :
 void BlockNonlinearForm::AddBdrFaceIntegrator(BlockNonlinearFormIntegrator *nfi,
                                               Array<int> &bdr_attr_marker)
 {
-   bfnfi.Append(nfi);
-   bfnfi_marker.Append(&bdr_attr_marker);
+   boundary_face_integs.Append(nfi);
+   boundary_face_integs_marker.Append(&bdr_attr_marker);
 }
 
 void BlockNonlinearForm::SetEssentialBC(
@@ -600,7 +769,7 @@ double BlockNonlinearForm::GetEnergyBlocked(const BlockVector &bx) const
       vdofs[i] = new Array<int>;
    }
 
-   if (dnfi.Size())
+   if (domain_integs.Size())
       for (int i = 0; i < fes[0]->GetNE(); ++i)
       {
          T = fes[0]->GetElementTransformation(i);
@@ -612,9 +781,9 @@ double BlockNonlinearForm::GetEnergyBlocked(const BlockVector &bx) const
             if (doftrans) {doftrans->InvTransformPrimal(*el_x[s]); }
          }
 
-         for (int k = 0; k < dnfi.Size(); ++k)
+         for (int k = 0; k < domain_integs.Size(); ++k)
          {
-            energy += dnfi[k]->GetElementEnergy(fe, *T, el_x_const);
+            energy += domain_integs[k]->GetElementEnergy(fe, *T, el_x_const);
          }
       }
 
@@ -625,12 +794,12 @@ double BlockNonlinearForm::GetEnergyBlocked(const BlockVector &bx) const
       delete vdofs[i];
    }
 
-   if (fnfi.Size())
+   if (interior_face_integs.Size())
    {
       MFEM_ABORT("TODO: add energy contribution from interior face terms");
    }
 
-   if (bfnfi.Size())
+   if (boundary_face_integs.Size())
    {
       MFEM_ABORT("TODO: add energy contribution from boundary face terms");
    }
@@ -668,7 +837,7 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
       vdofs2[s] = new Array<int>;
    }
 
-   if (dnfi.Size())
+   if (domain_integs.Size())
    {
       for (int i = 0; i < fes[0]->GetNE(); ++i)
       {
@@ -681,10 +850,10 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
             if (doftrans[s]) {doftrans[s]->InvTransformPrimal(*el_x[s]); }
          }
 
-         for (int k = 0; k < dnfi.Size(); ++k)
+         for (int k = 0; k < domain_integs.Size(); ++k)
          {
-            dnfi[k]->AssembleElementVector(fe, *T,
-                                           el_x_const, el_y);
+            domain_integs[k]->AssembleElementVector(fe, *T,
+                                                    el_x_const, el_y);
 
             for (int s=0; s<fes.Size(); ++s)
             {
@@ -696,7 +865,7 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
       }
    }
 
-   if (fnfi.Size())
+   if (interior_face_integs.Size())
    {
       Mesh *mesh = fes[0]->GetMesh();
       FaceElementTransformations *tr;
@@ -719,10 +888,10 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
                bx.GetBlock(s).GetSubVector(*(vdofs[s]), *el_x[s]);
             }
 
-            for (int k = 0; k < fnfi.Size(); ++k)
+            for (int k = 0; k < interior_face_integs.Size(); ++k)
             {
 
-               fnfi[k]->AssembleFaceVector(fe, fe2, *tr, el_x_const, el_y);
+               interior_face_integs[k]->AssembleFaceVector(fe, fe2, *tr, el_x_const, el_y);
 
                for (int s=0; s<fes.Size(); ++s)
                {
@@ -734,7 +903,7 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
       }
    }
 
-   if (bfnfi.Size())
+   if (boundary_face_integs.Size())
    {
       Mesh *mesh = fes[0]->GetMesh();
       FaceElementTransformations *tr;
@@ -742,14 +911,14 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
       Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
                                  mesh->bdr_attributes.Max() : 0);
       bdr_attr_marker = 0;
-      for (int k = 0; k < bfnfi.Size(); ++k)
+      for (int k = 0; k < boundary_face_integs.Size(); ++k)
       {
-         if (bfnfi_marker[k] == NULL)
+         if (boundary_face_integs_marker[k] == NULL)
          {
             bdr_attr_marker = 1;
             break;
          }
-         Array<int> &bdr_marker = *bfnfi_marker[k];
+         Array<int> &bdr_marker = *boundary_face_integs_marker[k];
          MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
                      "invalid boundary marker for boundary face integrator #"
                      << k << ", counting from zero");
@@ -776,12 +945,12 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
                bx.GetBlock(s).GetSubVector(*(vdofs[s]), *el_x[s]);
             }
 
-            for (int k = 0; k < bfnfi.Size(); ++k)
+            for (int k = 0; k < boundary_face_integs.Size(); ++k)
             {
-               if (bfnfi_marker[k] &&
-                   (*bfnfi_marker[k])[bdr_attr-1] == 0) { continue; }
+               if (boundary_face_integs_marker[k] &&
+                   (*boundary_face_integs_marker[k])[bdr_attr-1] == 0) { continue; }
 
-               bfnfi[k]->AssembleFaceVector(fe, fe2, *tr, el_x_const, el_y);
+               boundary_face_integs[k]->AssembleFaceVector(fe, fe2, *tr, el_x_const, el_y);
 
                for (int s=0; s<fes.Size(); ++s)
                {
@@ -886,7 +1055,7 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
       }
    }
 
-   if (dnfi.Size())
+   if (domain_integs.Size())
    {
       for (int i = 0; i < fes[0]->GetNE(); ++i)
       {
@@ -899,9 +1068,9 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
             if (doftrans[s]) {doftrans[s]->InvTransformPrimal(*el_x[s]); }
          }
 
-         for (int k = 0; k < dnfi.Size(); ++k)
+         for (int k = 0; k < domain_integs.Size(); ++k)
          {
-            dnfi[k]->AssembleElementGrad(fe, *T, el_x_const, elmats);
+            domain_integs[k]->AssembleElementGrad(fe, *T, el_x_const, elmats);
 
             for (int j=0; j<fes.Size(); ++j)
             {
@@ -920,7 +1089,7 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
       }
    }
 
-   if (fnfi.Size())
+   if (interior_face_integs.Size())
    {
       FaceElementTransformations *tr;
       Mesh *mesh = fes[0]->GetMesh();
@@ -941,9 +1110,9 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
             bx.GetBlock(s).GetSubVector(*vdofs[s], *el_x[s]);
          }
 
-         for (int k = 0; k < fnfi.Size(); ++k)
+         for (int k = 0; k < interior_face_integs.Size(); ++k)
          {
-            fnfi[k]->AssembleFaceGrad(fe, fe2, *tr, el_x_const, elmats);
+            interior_face_integs[k]->AssembleFaceGrad(fe, fe2, *tr, el_x_const, elmats);
             for (int j=0; j<fes.Size(); ++j)
             {
                for (int l=0; l<fes.Size(); ++l)
@@ -957,7 +1126,7 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
       }
    }
 
-   if (bfnfi.Size())
+   if (boundary_face_integs.Size())
    {
       FaceElementTransformations *tr;
       Mesh *mesh = fes[0]->GetMesh();
@@ -966,14 +1135,14 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
       Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
                                  mesh->bdr_attributes.Max() : 0);
       bdr_attr_marker = 0;
-      for (int k = 0; k < bfnfi.Size(); ++k)
+      for (int k = 0; k < boundary_face_integs.Size(); ++k)
       {
-         if (bfnfi_marker[k] == NULL)
+         if (boundary_face_integs_marker[k] == NULL)
          {
             bdr_attr_marker = 1;
             break;
          }
-         Array<int> &bdr_marker = *bfnfi_marker[k];
+         Array<int> &bdr_marker = *boundary_face_integs_marker[k];
          MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
                      "invalid boundary marker for boundary face integrator #"
                      << k << ", counting from zero");
@@ -1000,11 +1169,11 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
                bx.GetBlock(s).GetSubVector(*vdofs[s], *el_x[s]);
             }
 
-            for (int k = 0; k < bfnfi.Size(); ++k)
+            for (int k = 0; k < boundary_face_integs.Size(); ++k)
             {
-               if (bfnfi_marker[k] &&
-                   (*bfnfi_marker[k])[bdr_attr-1] == 0) { continue; }
-               bfnfi[k]->AssembleFaceGrad(fe, fe2, *tr, el_x_const, elmats);
+               if (boundary_face_integs_marker[k] &&
+                   (*boundary_face_integs_marker[k])[bdr_attr-1] == 0) { continue; }
+               boundary_face_integs[k]->AssembleFaceGrad(fe, fe2, *tr, el_x_const, elmats);
                for (int l=0; l<fes.Size(); ++l)
                {
                   for (int j=0; j<fes.Size(); ++j)
@@ -1109,19 +1278,19 @@ BlockNonlinearForm::~BlockNonlinearForm()
       delete ess_tdofs[i];
    }
 
-   for (int i = 0; i < dnfi.Size(); ++i)
+   for (int i = 0; i < domain_integs.Size(); ++i)
    {
-      delete dnfi[i];
+      delete domain_integs[i];
    }
 
-   for (int i = 0; i < fnfi.Size(); ++i)
+   for (int i = 0; i < interior_face_integs.Size(); ++i)
    {
-      delete fnfi[i];
+      delete interior_face_integs[i];
    }
 
-   for (int i = 0; i < bfnfi.Size(); ++i)
+   for (int i = 0; i < boundary_face_integs.Size(); ++i)
    {
-      delete bfnfi[i];
+      delete boundary_face_integs[i];
    }
 
 }

fem/nonlinearform.hpp

@@ -36,14 +36,19 @@ protected:
    FiniteElementSpace *fes; // not owned
 
    /// Set of Domain Integrators to be assembled (added).
-   Array<NonlinearFormIntegrator*> dnfi; // owned
+   Array<NonlinearFormIntegrator*> domain_integs; // owned
+   /// Element attribute marker (should be of length mesh->attributes)
+   /// Includes all by default.
+   /// 0 - ignore attribute
+   /// 1 - include attribute
+   Array<Array<int>*>             domain_integs_marker; // not owned
 
    /// Set of interior face Integrators to be assembled (added).
-   Array<NonlinearFormIntegrator*> fnfi; // owned
+   Array<NonlinearFormIntegrator*> interior_face_integs; // owned
 
    /// Set of boundary face Integrators to be assembled (added).
-   Array<NonlinearFormIntegrator*> bfnfi; // owned
-   Array<Array<int>*>              bfnfi_marker; // not owned
+   Array<NonlinearFormIntegrator*> boundary_face_integs; // owned
+   Array<Array<int>*>              boundary_face_integs_marker; // not owned
 
    mutable SparseMatrix *Grad, *cGrad; // owned
    /// Gradient Operator when not assembled as a matrix.
@@ -77,6 +82,20 @@ public:
         cP(dynamic_cast<const SparseMatrix*>(P))
    { }
 
+   /// Explicitly prohibit copy construction/assignment of NonlinearForms
+   NonlinearForm(const NonlinearForm &other) = delete;
+   NonlinearForm& operator=(const NonlinearForm &other) = delete;
+
+   /// Move constructor for NonlinearForm.
+   /** This constructor "steals" the owned data members from the @a other
+       NonlinearForm. */
+   NonlinearForm(NonlinearForm &&other);
+
+   /// Move assignment operator for NonlinearForm
+   /** This assignment first frees all owned data, then "steals" the owned data
+       members from the @a other NonlinearForm. */
+   NonlinearForm& operator=(NonlinearForm &&other);
+
    /// Set the desired assembly level. The default is AssemblyLevel::LEGACY.
    /** For nonlinear operators, the "matrix" assembly levels usually do not make
        sense, so only LEGACY, NONE (matrix-free) and PARTIAL are supported.
@@ -106,37 +125,44 @@ public:
    FiniteElementSpace *FESpace() { return fes; }
    const FiniteElementSpace *FESpace() const { return fes; }
 
-   /// Adds new Domain Integrator.
-   void AddDomainIntegrator(NonlinearFormIntegrator *nlfi)
-   { dnfi.Append(nlfi); }
+   /** @brief Adds new domain integrator.
+       @note Assumes ownership of @a nlfi. */
+   void AddDomainIntegrator(NonlinearFormIntegrator *nlfi);
+
+   /** @brief Adds new domain integrator restricted to certain elements
+       specified by the @a elem_marker.
+       @note Assumes ownership of @a nlfi. The array @a elem_marker is stored
+       internally as a pointer to the given Array<int> object. */
+   void AddDomainIntegrator(NonlinearFormIntegrator *nlfi,
+                            Array<int> &elem_marker);
 
    /// Access all integrators added with AddDomainIntegrator().
-   Array<NonlinearFormIntegrator*> *GetDNFI() { return &dnfi; }
-   const Array<NonlinearFormIntegrator*> *GetDNFI() const { return &dnfi; }
+   Array<NonlinearFormIntegrator*> *GetDNFI() { return &domain_integs; }
+   const Array<NonlinearFormIntegrator*> *GetDNFI() const { return &domain_integs; }
 
    /// Adds new Interior Face Integrator.
-   void AddInteriorFaceIntegrator(NonlinearFormIntegrator *nlfi)
-   { fnfi.Append(nlfi); }
+   void AddInteriorFaceIntegrator(NonlinearFormIntegrator *nlfi);
 
    /** @brief Access all interior face integrators added with
        AddInteriorFaceIntegrator(). */
    const Array<NonlinearFormIntegrator*> &GetInteriorFaceIntegrators() const
-   { return fnfi; }
+   { return interior_face_integs; }
 
-   /// Adds new Boundary Face Integrator.
-   void AddBdrFaceIntegrator(NonlinearFormIntegrator *nlfi)
-   { bfnfi.Append(nlfi); bfnfi_marker.Append(NULL); }
+   /** @brief Adds new boundary face integrator.
+       @note Assumes ownership of @a nlfi. */
+   void AddBdrFaceIntegrator(NonlinearFormIntegrator *nlfi);
 
-   /** @brief Adds new Boundary Face Integrator, restricted to specific boundary
-       attributes. */
-   void AddBdrFaceIntegrator(NonlinearFormIntegrator *nfi,
-                             Array<int> &bdr_marker)
-   { bfnfi.Append(nfi); bfnfi_marker.Append(&bdr_marker); }
+   /** @brief Adds new boundary face integrator restricted to specific boundary
+       attributes by the @a bdr_marker.
+       @note Assumes ownership of @a nlfi. The array @a bdr_marker is stored
+       internally as a pointer to the given Array<int> object. */
+   void AddBdrFaceIntegrator(NonlinearFormIntegrator *nlfi,
+                             Array<int> &bdr_marker);
 
    /** @brief Access all boundary face integrators added with
        AddBdrFaceIntegrator(). */
    const Array<NonlinearFormIntegrator*> &GetBdrFaceIntegrators() const
-   { return bfnfi; }
+   { return boundary_face_integs; }
 
    /// Specify essential boundary conditions.
    /** This method calls FiniteElementSpace::GetEssentialTrueDofs() and stores
@@ -226,14 +252,14 @@ protected:
    Array<FiniteElementSpace*> fes;
 
    /// Set of Domain Integrators to be assembled (added).
-   Array<BlockNonlinearFormIntegrator*> dnfi;
+   Array<BlockNonlinearFormIntegrator*> domain_integs;
 
    /// Set of interior face Integrators to be assembled (added).
-   Array<BlockNonlinearFormIntegrator*> fnfi;
+   Array<BlockNonlinearFormIntegrator*> interior_face_integs;
 
    /// Set of Boundary Face Integrators to be assembled (added).
-   Array<BlockNonlinearFormIntegrator*> bfnfi;
-   Array<Array<int>*>           bfnfi_marker;
+   Array<BlockNonlinearFormIntegrator*> boundary_face_integs;
+   Array<Array<int>*>           boundary_face_integs_marker;
 
    /** Auxiliary block-vectors for wrapping input and output vectors or holding
        GridFunction-like block-vector data (e.g. in parallel). */
@@ -298,15 +324,15 @@ public:
 
    /// Adds new Domain Integrator.
    void AddDomainIntegrator(BlockNonlinearFormIntegrator *nlfi)
-   { dnfi.Append(nlfi); }
+   { domain_integs.Append(nlfi); }
 
    /// Adds new Interior Face Integrator.
    void AddInteriorFaceIntegrator(BlockNonlinearFormIntegrator *nlfi)
-   { fnfi.Append(nlfi); }
+   { interior_face_integs.Append(nlfi); }
 
    /// Adds new Boundary Face Integrator.
    void AddBdrFaceIntegrator(BlockNonlinearFormIntegrator *nlfi)
-   { bfnfi.Append(nlfi); bfnfi_marker.Append(NULL); }
+   { boundary_face_integs.Append(nlfi); boundary_face_integs_marker.Append(NULL); }
 
    /** @brief Adds new Boundary Face Integrator, restricted to specific boundary
        attributes. */

fem/pbilinearform.cpp

@@ -19,6 +19,55 @@
 namespace mfem
 {
 
+ParBilinearForm::ParBilinearForm(ParBilinearForm &&other)
+   : BilinearForm(std::move(other)), pfes(other.pfes),
+     Xaux(other.pfes, other.Xaux.GetData()), Yaux(other.pfes, other.Yaux.GetData()),
+     Ytmp(std::move(other.Ytmp)), p_mat(other.p_mat), p_mat_e(other.p_mat_e),
+     keep_nbr_block(other.keep_nbr_block)
+{
+   other.Xaux.MakeRef(other.pfes, nullptr);
+   other.Yaux.MakeRef(other.pfes, nullptr);
+
+   p_mat.SetOperatorOwner();
+   other.p_mat.SetOperatorOwner(false);
+   other.p_mat.SetType(Operator::Hypre_ParCSR);
+
+   p_mat_e.SetOperatorOwner();
+   other.p_mat_e.SetOperatorOwner(false);
+   other.p_mat_e.SetType(Operator::Hypre_ParCSR);
+
+   other.keep_nbr_block = false;
+}
+
+ParBilinearForm& ParBilinearForm::operator=(ParBilinearForm &&other)
+{
+   if (this != &other)
+   {
+      BilinearForm::operator=(std::move(other));
+      pfes = other.pfes;
+      Xaux.MakeRef(other.pfes, other.Xaux.GetData());
+      other.Xaux.MakeRef(other.pfes, nullptr);
+      Yaux.MakeRef(other.pfes, other.Yaux.GetData());
+      other.Yaux.MakeRef(other.pfes, nullptr);
+
+      Ytmp = std::move(other.Ytmp);
+
+      p_mat = other.p_mat;
+      p_mat.SetOperatorOwner();
+      other.p_mat.SetOperatorOwner(false);
+      other.p_mat.SetType(Operator::Hypre_ParCSR);
+
+      p_mat_e = other.p_mat_e;
+      p_mat_e.SetOperatorOwner();
+      other.p_mat_e.SetOperatorOwner(false);
+      other.p_mat_e.SetType(Operator::Hypre_ParCSR);
+
+      keep_nbr_block = other.keep_nbr_block;
+      other.keep_nbr_block = false;
+   }
+   return *this;
+}
+
 void ParBilinearForm::pAllocMat()
 {
    int nbr_size = pfes->GetFaceNbrVSize();

fem/pbilinearform.hpp

@@ -44,13 +44,6 @@ protected:
 
    void AssembleSharedFaces(int skip_zeros = 1);
 
-private:
-   /// Copy construction is not supported; body is undefined.
-   ParBilinearForm(const ParBilinearForm &);
-
-   /// Copy assignment is not supported; body is undefined.
-   ParBilinearForm &operator=(const ParBilinearForm &);
-
 public:
    /// Creates parallel bilinear form associated with the FE space @a *pf.
    /** The pointer @a pf is not owned by the newly constructed object. */
@@ -71,6 +64,13 @@ public:
         p_mat(Operator::Hypre_ParCSR), p_mat_e(Operator::Hypre_ParCSR)
    { keep_nbr_block = false; }
 
+   /// Explicitly prohibit copy construction/assignment of ParBilinearForm
+   ParBilinearForm(const ParBilinearForm &) = delete;
+   ParBilinearForm &operator=(const ParBilinearForm &) = delete;
+
+   ParBilinearForm(ParBilinearForm &&other);
+   ParBilinearForm& operator=(ParBilinearForm &&other);
+
    /** When set to true and the ParBilinearForm has interior face integrators,
        the local SparseMatrix will include the rows (in addition to the columns)
        corresponding to face-neighbor dofs. The default behavior is to disregard

fem/plinearform.cpp

@@ -18,6 +18,17 @@
 namespace mfem
 {
 
+ParLinearForm& ParLinearForm::operator=(ParLinearForm &&other)
+{
+   if (this != &other)
+   {
+      LinearForm::operator=(std::move(other));
+      pfes = other.pfes;
+      other.pfes = nullptr;
+   }
+   return *this;
+}
+
 void ParLinearForm::Update(ParFiniteElementSpace *pf)
 {
    if (pf) { pfes = pf; }

fem/plinearform.hpp

@@ -28,10 +28,6 @@ class ParLinearForm : public LinearForm
 protected:
    ParFiniteElementSpace *pfes; ///< Points to the same object as #fes
 
-private:
-   /// Copy construction is not supported; body is undefined.
-   ParLinearForm(const ParLinearForm &);
-
 public:
    /** @brief Create an empty ParLinearForm without an associated
        ParFiniteElementSpace.
@@ -64,6 +60,9 @@ public:
    ParLinearForm(ParFiniteElementSpace *pf, ParLinearForm * plf)
       : LinearForm(pf, plf) { pfes = pf; }
 
+   /// Explicitly prohibit copy construction of ParLinearForm
+   ParLinearForm(const ParLinearForm &other) = delete;
+
    /// Copy assignment. Only the data of the base class Vector is copied.
    /** It is assumed that this object and @a rhs use ParFiniteElementSpace%s
        that have the same size.
@@ -73,6 +72,18 @@ public:
    ParLinearForm &operator=(const ParLinearForm &rhs)
    { return operator=((const Vector &)rhs); }
 
+   /// Move constructor for ParLinearForm.
+   /** This constructor "steals" the owned data members from the @a other
+       ParLinearForm. */
+   ParLinearForm(ParLinearForm &&other)
+      : LinearForm(std::move(other)), pfes(other.pfes)
+   { other.pfes = nullptr; }
+
+   /// Move assignment operator for ParLinearForm
+   /** This assignment first frees all owned data, then "steals" the owned data
+       members from the @a other ParLinearForm. */
+   ParLinearForm& operator=(ParLinearForm &&other);
+
    ParFiniteElementSpace *ParFESpace() const { return pfes; }
 
    /// Update the object according to the given new FE space @a *pf.

fem/pnonlinearform.cpp

@@ -27,13 +27,45 @@ ParNonlinearForm::ParNonlinearForm(ParFiniteElementSpace *pf)
    MFEM_VERIFY(!Serial(), "internal MFEM error");
 }
 
+ParNonlinearForm::ParNonlinearForm(ParNonlinearForm &&other)
+   : NonlinearForm(std::move(other)),
+     X((ParFiniteElementSpace *)other.fes, other.X.GetData()),
+     Y((ParFiniteElementSpace *)other.fes, other.Y.GetData()), pGrad(other.pGrad)
+{
+   other.X.MakeRef(other.fes, nullptr);
+   other.Y.MakeRef(other.fes, nullptr);
+
+   pGrad.SetOperatorOwner();
+   other.pGrad.SetOperatorOwner(false);
+   other.pGrad.SetType(Operator::Hypre_ParCSR);
+}
+
+ParNonlinearForm& ParNonlinearForm::operator=(ParNonlinearForm &&other)
+{
+   if (this != &other)
+   {
+      NonlinearForm::operator=(std::move(other));
+      X.MakeRef(other.fes, other.X.GetData());
+      other.X.MakeRef(other.fes, nullptr);
+      Y.MakeRef(other.fes, other.Y.GetData());
+      other.Y.MakeRef(other.fes, nullptr);
+
+      pGrad = other.pGrad;
+      pGrad.SetOperatorOwner();
+      other.pGrad.SetOperatorOwner(false);
+      other.pGrad.SetType(Operator::Hypre_ParCSR);
+
+   }
+   return *this;
+}
+
 double ParNonlinearForm::GetParGridFunctionEnergy(const Vector &x) const
 {
    double loc_energy, glob_energy;
 
    loc_energy = GetGridFunctionEnergy(x);
 
-   if (fnfi.Size())
+   if (interior_face_integs.Size())
    {
       MFEM_ABORT("TODO: add energy contribution from shared faces");
    }
@@ -48,7 +80,7 @@ void ParNonlinearForm::Mult(const Vector &x, Vector &y) const
 {
    NonlinearForm::Mult(x, y); // x --(P)--> aux1 --(A_local)--> aux2
 
-   if (fnfi.Size())
+   if (interior_face_integs.Size())
    {
       MFEM_VERIFY(!NonlinearForm::ext, "Not implemented (extensions + faces");
       // Terms over shared interior faces in parallel.
@@ -78,9 +110,9 @@ void ParNonlinearForm::Mult(const Vector &x, Vector &y) const
          X.GetSubVector(vdofs1, el_x.GetData());
          X.FaceNbrData().GetSubVector(vdofs2, el_x.GetData() + vdofs1.Size());
 
-         for (int k = 0; k < fnfi.Size(); k++)
+         for (int k = 0; k < interior_face_integs.Size(); k++)
          {
-            fnfi[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
+            interior_face_integs[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
             aux2.AddElementVector(vdofs1, el_y.GetData());
          }
       }
@@ -116,7 +148,7 @@ Operator &ParNonlinearForm::GetGradient(const Vector &x) const
 
    OperatorHandle dA(pGrad.Type()), Ph(pGrad.Type());
 
-   if (fnfi.Size() == 0)
+   if (interior_face_integs.Size() == 0)
    {
       dA.MakeSquareBlockDiag(pfes->GetComm(), pfes->GlobalVSize(),
                              pfes->GetDofOffsets(), Grad);
@@ -252,7 +284,7 @@ void ParBlockNonlinearForm::Mult(const Vector &x, Vector &y) const
 
    BlockNonlinearForm::MultBlocked(xs, ys);
 
-   if (fnfi.Size() > 0)
+   if (interior_face_integs.Size() > 0)
    {
       MFEM_ABORT("TODO: assemble contributions from shared face terms");
    }
@@ -332,7 +364,7 @@ BlockOperator & ParBlockNonlinearForm::GetGradient(const Vector &x) const
 
    GetLocalGradient(x); // gradients are stored in 'Grads'
 
-   if (fnfi.Size() > 0)
+   if (interior_face_integs.Size() > 0)
    {
       MFEM_ABORT("TODO: assemble contributions from shared face terms");
    }

fem/pnonlinearform.hpp

@@ -32,6 +32,13 @@ protected:
 public:
    ParNonlinearForm(ParFiniteElementSpace *pf);
 
+   /// Explicitly prohibit copy construction/assignment of ParNonlinearForm
+   ParNonlinearForm(const ParNonlinearForm &other) = delete;
+   ParNonlinearForm& operator=(const ParNonlinearForm &other) = delete;
+
+   ParNonlinearForm(ParNonlinearForm &&other);
+   ParNonlinearForm& operator=(ParNonlinearForm &&other);
+
    ParFiniteElementSpace *ParFESpace() const
    { return (ParFiniteElementSpace *)fes; }