Merge remote-tracking branch 'upstream/feature/partial_ess_bcs' into exaconstit-csr

Removed the more or less duplicated FESpace::GetEssentialVDofs and moved the necessary book-keeping portions of things to GetEssentialTrueDofs
Additionally, I added an overwrite optional variable (defaulted to true) to GetEssentialVDofs in-order to allow one to overwrite the values of ess_vdofs if they'd like.
I still need to add some unit tests for the functionality though.

data/cube-hex-ro.mesh

@@ -0,0 +1,86 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+#
+
+dimension
+3
+
+elements
+8
+1 5  0   1   10  9   3   4   13  12 
+1 5  1   2   11  10  4   5   14  13 
+1 5  3   4   13  12  6   7   16  15
+1 5  4   5   14  13  7   8   17  16
+1 5  9   10  19  18  12  13  22  21
+1 5  10  11  20  19  13  14  23  22
+1 5  12  13  22  21  15  16  25  24
+1 5  13  14  23  22  16  17  26  25
+
+#
+
+boundary
+24
+1 3 1  0  9  10 
+1 3 2  1  10 11 
+1 3 10 9  18 19 
+1 3 11 10 19 20 
+2 3 0  3  12 9
+2 3 9  12 21 18
+2 3 3  6  15 12
+2 3 12 15 24 21
+3 3 0  1  4  3
+3 3 1  2  5  4
+3 3 3  4  7  6
+3 3 4  5  8  7 
+4 3 6  7  16 15
+4 3 7  8  17 16
+4 3 15 16 25 24
+4 3 16 17 26 25 
+5 3 18 21 22 19 
+5 3 19 22 23 20 
+5 3 21 24 25 22 
+5 3 22 25 26 23 
+6 3 2  11 14 5
+6 3 11 20 23 14
+6 3 5  14 17 8 
+6 3 14 23 26 17
+
+vertices
+27
+3
+0.0 0.0 0.0
+0.5 0.0 0.0
+1.0 0.0 0.0
+0.0 0.0 0.5
+0.5 0.0 0.5
+1.0 0.0 0.5
+0.0 0.0 1.0
+0.5 0.0 1.0
+1.0 0.0 1.0
+0.0 0.5 0.0
+0.5 0.5 0.0
+1.0 0.5 0.0
+0.0 0.5 0.5
+0.5 0.5 0.5
+1.0 0.5 0.5
+0.0 0.5 1.0
+0.5 0.5 1.0
+1.0 0.5 1.0
+0.0 1.0 0.0
+0.5 1.0 0.0
+1.0 1.0 0.0
+0.0 1.0 0.5
+0.5 1.0 0.5
+1.0 1.0 0.5
+0.0 1.0 1.0
+0.5 1.0 1.0
+1.0 1.0 1.0

data/cube-hex.mesh

@@ -0,0 +1,84 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+#
+
+dimension
+3
+
+elements
+8
+1 5 0  1  4  3  9  10 13 12 
+1 5 1  2  5  4  10 11 14 13
+1 5 9  10 13 12 18 19 22 21
+1 5 10 11 14 13 19 20 23 22
+1 5 3  4  7  6  12 13 16 15
+1 5 4  5  8  7  13 14 17 16
+1 5 12 13 16 15 21 22 25 24
+1 5 13 14 17 16 22 23 26 25
+
+boundary
+24
+1 3 0  1  10 9 
+1 3 1  2  11 10 
+1 3 9  10 19 18 
+1 3 10 11 20 19 
+3 3 0  3  4  1 
+3 3 1  4  5  2 
+3 3 3  6  7  4 
+3 3 4  7  8  5 
+3 3 18 19 22 21 
+3 3 19 20 23 22 
+3 3 21 22 25 24 
+3 3 22 23 26 25 
+3 3 2  5  14 11 
+3 3 11 14 23 20 
+3 3 5  8  17 14 
+3 3 14 17 26 23 
+3 3 0  9  12 3 
+3 3 9  18 21 12 
+3 3 3  12 15 6 
+3 3 12 21 24 15 
+2 3 6  15 16 7 
+2 3 7  16 17 8 
+2 3 15 24 25 16 
+2 3 16 25 26 17 
+
+vertices
+27
+3
+0.0 0.0 0.0
+0.5 0.0 0.0
+1.0 0.0 0.0
+0.0 0.0 0.5
+0.5 0.0 0.5
+1.0 0.0 0.5
+0.0 0.0 1.0
+0.5 0.0 1.0
+1.0 0.0 1.0
+0.0 0.5 0.0
+0.5 0.5 0.0
+1.0 0.5 0.0
+0.0 0.5 0.5
+0.5 0.5 0.5
+1.0 0.5 0.5
+0.0 0.5 1.0
+0.5 0.5 1.0
+1.0 0.5 1.0
+0.0 1.0 0.0
+0.5 1.0 0.0
+1.0 1.0 0.0
+0.0 1.0 0.5
+0.5 1.0 0.5
+1.0 1.0 0.5
+0.0 1.0 1.0
+0.5 1.0 1.0
+1.0 1.0 1.0

fem/coefficient.cpp

@@ -517,6 +517,30 @@ void VectorRestrictedCoefficient::Eval(
    }
 }
 
+void VectorFunctionRestrictedCoefficient::Eval(Vector &V,
+                                               ElementTransformation &T,
+                                               const IntegrationPoint &ip)
+{
+   double x[3];
+   Vector transip(x, 3);
+
+   T.Transform(ip, transip);
+
+   V.SetSize(vdim);
+   if (active_attr[T.Attribute-1])
+   {
+      (*TDFunction)(T.Attribute, V);
+   }
+   else
+   {
+      V = 0.0;
+   }
+   if (Q)
+   {
+      V *= Q->Eval(T, ip, GetTime());
+   }
+}
+
 void PWMatrixCoefficient::InitMap(const Array<int> & attr,
                                   const Array<MatrixCoefficient*> & coefs)
 {

fem/coefficient.hpp

@@ -439,12 +439,19 @@ public:
    /// Initialize the VectorCoefficient with vector dimension @a vd.
    VectorCoefficient(int vd) { vdim = vd; time = 0.; }
 
+   // add constructor with no argument list (this is for the
+   // QuadratureVectorFunctionCoefficient extended class), srw
+   VectorCoefficient() { vdim = 0; time = 0.; }
+
    /// Set the time for time dependent coefficients
    virtual void SetTime(double t) { time = t; }
 
    /// Get the time for time dependent coefficients
    double GetTime() { return time; }
 
+   /// Set the dimension of the vector.
+   void SetVDim(int dim) { vdim = dim; }
+
    /// Returns dimension of the vector.
    int GetVDim() { return vdim; }
 
@@ -870,6 +877,37 @@ public:
                      const IntegrationRule &ir);
 };
 
+/// VectorFunctionCoefficient defined on a subset of domain or boundary attributes, srw
+class VectorFunctionRestrictedCoefficient : public VectorCoefficient
+{
+private:
+   void (*TDFunction)(int, Vector &);
+   const Array<int> &active_attr;
+   const Array2D<double> &active_scale;
+   Coefficient *Q;
+
+public:
+   /// Construct a time-dependent vector coefficient from a C-function
+   VectorFunctionRestrictedCoefficient(int dim,
+                                       void (*TDF)(int, Vector &),
+                                       Array<int> &attr,
+                                       Array2D<double> &scale,
+                                       Coefficient *q = NULL)
+      : VectorCoefficient(dim), Q(q), active_attr(attr), active_scale(scale)
+   {
+      TDFunction = TDF;
+   }
+
+   using VectorCoefficient::Eval;
+   virtual void Eval(Vector &V, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+
+   virtual ~VectorFunctionRestrictedCoefficient() { }
+
+   const Array<int> &GetActiveAttr() { return active_attr; }
+   const Array2D<double> &GetActiveScale() { return active_scale; }
+
+};
 typedef VectorCoefficient DiagonalMatrixCoefficient;
 
 /// Base class for Matrix Coefficients that optionally depend on time and space.

fem/fespace.cpp

@@ -494,12 +494,21 @@ static void mark_dofs(const Array<int> &dofs, Array<int> &mark_array)
 
 void FiniteElementSpace::GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
                                            Array<int> &ess_vdofs,
-                                           int component) const
+                                           int component,
+                                           bool overwrite) const
 {
    Array<int> vdofs, dofs;
 
-   ess_vdofs.SetSize(GetVSize());
-   ess_vdofs = 0;
+   if (overwrite)
+   {
+      ess_vdofs.SetSize(GetVSize());
+      ess_vdofs = 0;
+   }
+   else
+   {
+      MFEM_ASSERT(ess_vdofs.Size() == GetVSize(),
+                  "ess_vdofs size is not equal to FESpaces GetVSize().");
+   }
 
    for (int i = 0; i < GetNBE(); i++)
    {
@@ -579,6 +588,47 @@ void FiniteElementSpace::GetEssentialTrueDofs(const Array<int> &bdr_attr_is_ess,
    MarkerToList(ess_tdofs, ess_tdof_list);
 }
 
+void FiniteElementSpace::GetEssentialTrueDofs(const Array<int> &bdr_attr_is_ess,
+                                              Array<int> &ess_tdof_list,
+                                              const Array2D<bool> &component)
+{
+   MFEM_ASSERT(component.NumCols() == vdim,
+               "Number of columns of component was not equal to FESpace vdim");
+   MFEM_ASSERT(component.NumRows() == bdr_attr_is_ess.Size(),
+               "Number of rows of component was not equal to bdr_attr_is_ess.Size()");
+
+   Array<int> ess_vdofs, ess_tdofs, bdr_attr_is_ess_single_comp;
+   bdr_attr_is_ess_single_comp.SetSize(bdr_attr_is_ess.Size());
+
+   for (int i = 0; i < vdim; i++)
+   {
+      // Only overwrite ess_vdofs on first iteration
+      // all other iterations we want to preserve values of
+      // ess_vdofs.
+      const bool overwrite = (i == 0) ? true : false;
+      bdr_attr_is_ess_single_comp = 0;
+      for (int j = 0; j < bdr_attr_is_ess.Size(); j++)
+      {
+         if (bdr_attr_is_ess[j] && component(j, i))
+         {
+            bdr_attr_is_ess_single_comp[j] = bdr_attr_is_ess[j];
+         }
+      }
+      GetEssentialVDofs(bdr_attr_is_ess_single_comp, ess_vdofs, i, overwrite);
+   }
+
+   const SparseMatrix *R = GetConformingRestriction();
+   if (!R)
+   {
+      ess_tdofs.MakeRef(ess_vdofs);
+   }
+   else
+   {
+      R->BooleanMult(ess_vdofs, ess_tdofs);
+   }
+   MarkerToList(ess_tdofs, ess_tdof_list);
+}
+
 void FiniteElementSpace::GetBoundaryTrueDofs(Array<int> &boundary_dofs,
                                              int component)
 {

fem/fespace.hpp

@@ -775,11 +775,16 @@ public:
 
    /** @brief Mark degrees of freedom associated with boundary elements with
        the specified boundary attributes (marked in 'bdr_attr_is_ess').
+
        For spaces with 'vdim' > 1, the 'component' parameter can be used
-       to restricts the marked vDOFs to the specified component. */
+       to restricts the marked vDOFs to the specified component.
+       If overwrite is set to false then values in ess_vdofs are preserved
+       and not reset. However, the assumption here is that ess_vdofs is set to
+       the correct size already.*/
    virtual void GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
                                   Array<int> &ess_vdofs,
-                                  int component = -1) const;
+                                  int component = -1,
+                                  bool overwrite = true) const;
 
    /** @brief Get a list of essential true dofs, ess_tdof_list, corresponding to the
        boundary attributes marked in the array bdr_attr_is_ess.
@@ -789,6 +794,20 @@ public:
                                      Array<int> &ess_tdof_list,
                                      int component = -1);
 
+
+   /** @brief Get a list of essential true dofs, ess_tdof_list, corresponding to the
+       boundary attributes marked in the array bdr_attr_is_ess.
+
+       For spaces with 'vdim' > 1, the 'component' array can be used
+       to restricts the marked tDOFs per boundary to the specified components.
+       If vdim > 1 then one can specify per boundary attribute which components
+       on a boundary are essential by assigning a value of true to its location
+       in the component array.
+       The component has dimensions number of boundary attributes x vdim. */
+   virtual void GetEssentialTrueDofs(const Array<int> &bdr_attr_is_ess,
+                                     Array<int> &ess_tdof_list,
+                                     const Array2D<bool> &component);
+
    /** @brief Get a list of all boundary true dofs, @a boundary_dofs. For spaces
        with 'vdim' > 1, the 'component' parameter can be used to restricts the
        marked tDOFs to the specified component. Equivalent to

fem/gridfunc.cpp

@@ -2025,6 +2025,7 @@ void GridFunction::AccumulateAndCountZones(Coefficient &coeff,
    // Local interpolation
    Array<int> vdofs;
    Vector vals;
+   this->HostReadWrite();
    *this = 0.0;
 
    HostReadWrite();
@@ -2066,6 +2067,7 @@ void GridFunction::AccumulateAndCountZones(VectorCoefficient &vcoeff,
    // Local interpolation
    Array<int> vdofs;
    Vector vals;
+   this->HostReadWrite();
    *this = 0.0;
 
    HostReadWrite();
@@ -2672,6 +2674,100 @@ void GridFunction::ProjectBdrCoefficient(Coefficient *coeff[], Array<int> &attr)
 #endif
 }
 
+void GridFunction::ProjectBdrCoefficient(VectorFunctionCoefficient &vfcoeff,
+                                         Array<int> &attr)
+{
+   int i, j, fdof, d, ind, vdim;
+   Vector val;
+   const FiniteElement *fe;
+   ElementTransformation *transf;
+   Array<int> vdofs;
+
+   this->HostReadWrite();
+
+   vdim = fes->GetVDim();
+   // loop over boundary elements
+   for (i = 0; i < fes->GetNBE(); i++)
+   {
+      // if boundary attribute is 1 (Dirichlet)
+      if (attr[fes->GetBdrAttribute(i) - 1])
+      {
+         fe = fes->GetBE(i);
+         fdof = fe->GetDof();
+         transf = fes->GetBdrElementTransformation(i);
+         const IntegrationRule &ir = fe->GetNodes();
+         fes->GetBdrElementVDofs(i, vdofs);
+
+         // loop over dofs
+         for (j = 0; j < fdof; j++)
+         {
+            const IntegrationPoint &ip = ir.IntPoint(j);
+            transf->SetIntPoint(&ip);
+
+            vfcoeff.Eval(val, *transf, ip);
+
+            // loop over vector dimensions
+            for (d = 0; d < vdim; d++)
+            {
+               ind = vdofs[fdof*d+j];
+               if ( (ind = vdofs[fdof*d+j]) < 0 )
+               {
+                  val(d) = -val(d), ind = -1-ind;
+               }
+               (*this)(ind) = val(d);
+            }
+         }
+      }
+   }
+
+}
+
+void GridFunction::ProjectBdrCoefficient(VectorFunctionRestrictedCoefficient &vfcoeff)
+
+{
+   int i, j, fdof, d, ind, vdim;
+   Vector val;
+   const FiniteElement *fe;
+   ElementTransformation *transf;
+   Array<int> vdofs;
+   const Array<int> &active_attr = vfcoeff.GetActiveAttr();
+   const Array2D<double> &scale = vfcoeff.GetActiveScale();
+
+   this->HostReadWrite();
+   vdim = fes->GetVDim();
+   // loop over boundary elements
+   for (i = 0; i < fes->GetNBE(); i++) {
+      // if boundary attribute is 1 (Dirichlet)
+      if (active_attr[fes->GetBdrAttribute(i) - 1]) {
+         // instantiate a BC object
+         fe = fes->GetBE(i);
+         fdof = fe->GetDof();
+         transf = fes->GetBdrElementTransformation(i);
+         const IntegrationRule &ir = fe->GetNodes();
+         fes->GetBdrElementVDofs(i, vdofs);
+
+         // loop over dofs
+         for (j = 0; j < fdof; j++) {
+            const IntegrationPoint &ip = ir.IntPoint(j);
+            transf->SetIntPoint(&ip);
+            vfcoeff.Eval(val, *transf, ip);
+            // loop over vector dimensions
+            for (d = 0; d < vdim; d++) {
+               // check if the vector component (i.e. dof) is not constrained by a
+               // partial essential BC
+               if (scale(transf->Attribute - 1, d) > 0.0) {
+                  ind = vdofs[fdof * d + j];
+                  if ( (ind = vdofs[fdof * d + j]) < 0) {
+                     val(d) = -val(d), ind = -1 - ind;
+                  }
+                  (*this)(ind) = val(d); // placing computed value in grid function
+               }
+            }
+         }
+      }
+   }
+}
+
 void GridFunction::ProjectBdrCoefficientNormal(
    VectorCoefficient &vcoeff, Array<int> &bdr_attr)
 {

fem/gridfunc.hpp

@@ -463,6 +463,20 @@ public:
        will not be touched. */
    virtual void ProjectBdrCoefficient(Coefficient *coeff[], Array<int> &attr);
 
+   /** add ProjectBdrCoefficient() taking a VectorFunctionCoefficient and projecting
+       onto the grid function, srw. */
+   void ProjectBdrCoefficient(VectorFunctionCoefficient &vfcoeff,
+                              Array<int> &attr);
+
+   /** add ProjectBdrCoefficient() taking a VectorFunctionRestrictedCoefficient, srw
+       This is the same as the function above taking a VectorFunctionCoefficient, it
+       is just that the function coefficient is now restricted. This function
+       specifically projects values associated with ExaConstits essential BC
+       management and input. This routine does not, at this time, project a
+       VectorFunctionRestrictedCoefficient onto all grid function dofs.
+       This function should be renamed to reflect this. */
+   void ProjectBdrCoefficient(VectorFunctionRestrictedCoefficient &vfcoeff);
+
    /** Project the normal component of the given VectorCoefficient on
        the boundary. Only boundary attributes that are marked in
        'bdr_attr' are projected. Assumes RT-type VectorFE GridFunction. */
@@ -724,7 +738,6 @@ public:
    virtual ~GridFunction() { Destroy(); }
 };
 
-
 /** Overload operator<< for std::ostream and GridFunction; valid also for the
     derived class ParGridFunction */
 std::ostream &operator<<(std::ostream &out, const GridFunction &sol);
@@ -928,7 +941,6 @@ public:
 /// Overload operator<< for std::ostream and QuadratureFunction.
 std::ostream &operator<<(std::ostream &out, const QuadratureFunction &qf);
 
-
 double ZZErrorEstimator(BilinearFormIntegrator &blfi,
                         GridFunction &u,
                         GridFunction &flux,
@@ -1042,6 +1054,7 @@ inline void QuadratureFunction::SetSpace(QuadratureSpace *qspace_,
 
 inline void QuadratureFunction::GetElementValues(int idx, Vector &values)
 {
+   // element offset is the number of quadrature points for that element
    const int s_offset = qspace->element_offsets[idx];
    const int sl_size = qspace->element_offsets[idx+1] - s_offset;
    values.NewDataAndSize(data + vdim*s_offset, vdim*sl_size);

fem/nonlinearform.cpp

@@ -53,6 +53,22 @@ void NonlinearForm::SetEssentialBC(const Array<int> &bdr_attr_is_ess,
    }
 }
 
+void NonlinearForm::SetEssentialBC(const Array<int> &bdr_attr_is_ess,
+                                   const Array2D<bool> &bdr_component,
+                                   Vector *rhs)
+{
+   // virtual call, works in parallel too
+   fes->GetEssentialTrueDofs(bdr_attr_is_ess, ess_tdof_list, bdr_component);
+
+   if (rhs)
+   {
+      for (int i = 0; i < ess_tdof_list.Size(); i++)
+      {
+         (*rhs)(ess_tdof_list[i]) = 0.0;
+      }
+   }
+}
+
 void NonlinearForm::SetEssentialVDofs(const Array<int> &ess_vdofs_list)
 {
    if (!P)

fem/nonlinearform.hpp

@@ -146,8 +146,26 @@ public:
        have zero entries at the essential true dofs. */
    void SetEssentialBC(const Array<int> &bdr_attr_is_ess, Vector *rhs = NULL);
 
+   /** @brief Specify essential boundary conditions.
+
+       For spaces with 'vdim' > 1, the 'bdr_component' array can be used
+       to restricts the marked tDOFs per boundary to the specified components.
+       If vdim > 1 then one can specify per boundary attribute which components
+       on a boundary are essential by assigning a value of true to its location
+       in the bdr_component array.
+       The bdr_component has dimensions number of boundary attributes x vdim
+
+       This method calls FiniteElementSpace::GetEssentialTrueDofs() and stores
+       the result internally for use by other methods. If the @a rhs pointer is
+       not NULL, its essential true dofs will be set to zero. This makes it
+       "compatible" with the output vectors from the Mult() method which also
+       have zero entries at the essential true dofs. */
+   void SetEssentialBC(const Array<int> &bdr_attr_is_ess,
+                       const Array2D<bool> &bdr_component,
+                       Vector *rhs);
+
    /// Specify essential boundary conditions.
-   /** Use either SetEssentialBC() or SetEssentialTrueDofs() if possible. */
+   /** @deprecated Use either SetEssentialBC() or SetEssentialTrueDofs(). */
    void SetEssentialVDofs(const Array<int> &ess_vdofs_list);
 
    /// Specify essential boundary conditions.
@@ -179,6 +197,7 @@ public:
        vectors, i.e. their size must be fes->GetTrueVSize(). */
    virtual void Mult(const Vector &x, Vector &y) const;
 
+
    /** @brief Compute the gradient Operator of the NonlinearForm corresponding
        to the state @a x. */
    /** Any previously specified essential boundary conditions will be

fem/nonlininteg.cpp

@@ -80,6 +80,15 @@ void NonlinearFormIntegrator::AssembleElementVector(
               " is not overloaded!");
 }
 
+void NonlinearFormIntegrator::AssembleElementVector(
+   const FiniteElement &el, ElementTransformation &Ttr_beg,
+   ElementTransformation &Ttr_end,
+   const Vector &elfun, Vector &elvect, const Vector &elvel)
+{
+   mfem_error("NonlinearFormIntegrator::AssembleElementVector"
+              " is not overloaded!");
+}
+
 void NonlinearFormIntegrator::AssembleFaceVector(
    const FiniteElement &el1, const FiniteElement &el2,
    FaceElementTransformations &Tr, const Vector &elfun, Vector &elvect)
@@ -105,6 +114,20 @@ void NonlinearFormIntegrator::AssembleFaceGrad(
               " is not overloaded!");
 }
 
+void NonlinearFormIntegrator::AssembleGradPA(
+   const FiniteElementSpace &fes)
+{
+   mfem_error("NonlinearFormIntegrator::AssembleGradPA"
+              " is not overloaded!");
+}
+
+void NonlinearFormIntegrator::AssembleEA(
+   const FiniteElementSpace &fes, Vector &emat)
+{
+   mfem_error("NonlinearFormIntegrator::AssembleEA"
+              " is not overloaded!");
+}
+
 double NonlinearFormIntegrator::GetElementEnergy(
    const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun)
 {
@@ -166,7 +189,6 @@ double BlockNonlinearFormIntegrator::GetElementEnergy(
    return 0.0;
 }
 
-
 double InverseHarmonicModel::EvalW(const DenseMatrix &J) const
 {
    Z.SetSize(J.Width());

fem/nonlininteg.hpp

@@ -57,6 +57,16 @@ public:
                                       ElementTransformation &Tr,
                                       const Vector &elfun, Vector &elvect);
 
+   //This function really only be used whenever UMATs are being used,
+   //since UMATs require a velocity based solution vector,
+   //However, they depend on a lot of variables based on displacements
+   //to evolve their models. So, we need to be able to take
+   //derivatives using the beginning and the end time step mesh nodes.
+   //The other models should be using the one above this
+   virtual void AssembleElementVector(const FiniteElement &el,
+                                      ElementTransformation &Ttr_beg,
+                                      ElementTransformation &Ttr_end,
+                                      const Vector &elfun, Vector &elvect, const Vector &elvel);
    /// @brief Perform the local action of the NonlinearFormIntegrator resulting
    /// from a face integral term.
    virtual void AssembleFaceVector(const FiniteElement &el1,
@@ -112,6 +122,14 @@ public:
        This method can be called only after the method AssemblePA() has been
        called. */
    virtual void AddMultPA(const Vector &x, Vector &y) const;
+   
+   /** @brief Prepare the integrator for partial assembly (PA) gradient
+       evaluations on the given FE space @a fes */
+   /** The result of the partial assembly is stored internally so that it can be
+       used later in the methods AddMultGradPA() and AssembleGradDiagonalPA().
+       The state Vector @a x is an E-vector. */                                                           
+   virtual void AssembleGradPA(const FiniteElementSpace &fes);
+   virtual void AssembleEA(const FiniteElementSpace &fes, Vector &emat);
 
    /// Method for partially assembled gradient action.
    /** All arguments are E-vectors. This method can be called only after the
@@ -188,17 +206,16 @@ public:
    virtual ~BlockNonlinearFormIntegrator() { }
 };
 
-
-/// Abstract class for hyperelastic models
-class HyperelasticModel
+/// Abstract class for nonlinear models
+class NonlinearModel
 {
 protected:
    ElementTransformation *Ttr; /**< Reference-element to target-element
                                     transformation. */
 
 public:
-   HyperelasticModel() : Ttr(NULL) { }
-   virtual ~HyperelasticModel() { }
+   NonlinearModel() : Ttr(NULL) { }
+   virtual ~NonlinearModel() { }
 
    /// A reference-element to target-element transformation that can be used to
    /// evaluate Coefficient%s.
@@ -206,11 +223,6 @@ public:
        point of interest. */
    void SetTransformation(ElementTransformation &Ttr_) { Ttr = &Ttr_; }
 
-   /** @brief Evaluate the strain energy density function, W = W(Jpt).
-       @param[in] Jpt  Represents the target->physical transformation
-                       Jacobian matrix. */
-   virtual double EvalW(const DenseMatrix &Jpt) const = 0;
-
    /** @brief Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
        @param[in] Jpt  Represents the target->physical transformation
                        Jacobian matrix.
@@ -235,6 +247,18 @@ public:
 };
 
 
+/// Abstract class for hyperelastic models
+class HyperelasticModel : public NonlinearModel
+{
+public:
+
+   /** @brief Evaluate the strain energy density function, W = W(Jpt).
+       @param[in] Jpt  Represents the target->physical transformation
+                       Jacobian matrix. */
+   virtual double EvalW(const DenseMatrix &Jpt) const = 0;
+
+};
+
 /** Inverse-harmonic hyperelastic model with a strain energy density function
     given by the formula: W(J) = (1/2) det(J) Tr((J J^t)^{-1}) where J is the
     deformation gradient. */

fem/pfespace.cpp

@@ -1020,9 +1020,11 @@ void ParFiniteElementSpace::Synchronize(Array<int> &ldof_marker) const
 
 void ParFiniteElementSpace::GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
                                               Array<int> &ess_dofs,
-                                              int component) const
+                                              int component,
+                                              bool overwrite) const
 {
-   FiniteElementSpace::GetEssentialVDofs(bdr_attr_is_ess, ess_dofs, component);
+   FiniteElementSpace::GetEssentialVDofs(bdr_attr_is_ess, ess_dofs, component,
+                                         overwrite);
 
    // Make sure that processors without boundary elements mark
    // their boundary dofs (if they have any).
@@ -1059,6 +1061,58 @@ void ParFiniteElementSpace::GetEssentialTrueDofs(const Array<int>
    MarkerToList(true_ess_dofs, ess_tdof_list);
 }
 
+void ParFiniteElementSpace::GetEssentialTrueDofs(const Array<int>
+                                                 &bdr_attr_is_ess,
+                                                 Array<int> &ess_tdof_list,
+                                                 const Array2D<bool> &component)
+{
+   MFEM_ASSERT(component.NumCols() == vdim,
+               "Number of columns of component was not equal to ParFESpace vdim");
+   MFEM_ASSERT(component.NumRows() == bdr_attr_is_ess.Size(),
+               "Number of rows of component was not equal to bdr_attr_is_ess.Size()");
+
+   Array<int> ess_dofs, true_ess_dofs, bdr_attr_is_ess_single_comp;
+   bdr_attr_is_ess_single_comp.SetSize(bdr_attr_is_ess.Size());
+
+   for (int i = 0; i < vdim; i++)
+   {
+      // Only overwrite ess_vdofs on first iteration
+      // all other iterations we want to preserve values of
+      // ess_vdofs.
+      const bool overwrite = (i == 0) ? true : false;
+      bdr_attr_is_ess_single_comp = 0;
+      for (int j = 0; j < bdr_attr_is_ess.Size(); j++)
+      {
+         if (bdr_attr_is_ess[j] && component(j, i))
+         {
+            bdr_attr_is_ess_single_comp[j] = bdr_attr_is_ess[j];
+         }
+      }
+      GetEssentialVDofs(bdr_attr_is_ess_single_comp, ess_dofs, i, overwrite);
+   }
+
+   GetRestrictionMatrix()->BooleanMult(ess_dofs, true_ess_dofs);
+
+#ifdef MFEM_DEBUG
+   // Verify that in boolean arithmetic: P^T ess_dofs = R ess_dofs.
+   Array<int> true_ess_dofs2(true_ess_dofs.Size());
+   HypreParMatrix *Pt = Dof_TrueDof_Matrix()->Transpose();
+   const int *ess_dofs_data = ess_dofs.HostRead();
+   Pt->BooleanMult(1, ess_dofs_data, 0, true_ess_dofs2);
+   delete Pt;
+   int counter = 0;
+   const int *ted = true_ess_dofs.HostRead();
+   for (int i = 0; i < true_ess_dofs.Size(); i++)
+   {
+      if (bool(ted[i]) != bool(true_ess_dofs2[i])) { counter++; }
+   }
+   MFEM_VERIFY(counter == 0, "internal MFEM error: counter = " << counter
+               << ", rank = " << MyRank);
+#endif
+
+   MarkerToList(true_ess_dofs, ess_tdof_list);
+}
+
 int ParFiniteElementSpace::GetLocalTDofNumber(int ldof) const
 {
    if (Nonconforming())

fem/pfespace.hpp

@@ -350,10 +350,18 @@ public:
        "partially conforming") space. */
    void Synchronize(Array<int> &ldof_marker) const;
 
-   /// Determine the boundary degrees of freedom
+   /** @brief Mark degrees of freedom associated with boundary elements with
+       the specified boundary attributes (marked in 'bdr_attr_is_ess').
+
+       For spaces with 'vdim' > 1, the 'component' parameter can be used
+       to restricts the marked vDOFs to the specified component.
+       If overwrite is set to false then values in ess_vdofs are preserved
+       and not reset. However, the assumption here is that ess_vdofs is set to
+       the correct size already.*/
    virtual void GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
                                   Array<int> &ess_dofs,
-                                  int component = -1) const;
+                                  int component = -1,
+                                  bool overwrite = true) const;
 
    /** Get a list of essential true dofs, ess_tdof_list, corresponding to the
        boundary attributes marked in the array bdr_attr_is_ess. */
@@ -361,6 +369,19 @@ public:
                                      Array<int> &ess_tdof_list,
                                      int component = -1);
 
+   /** @brief Get a list of essential true dofs, ess_tdof_list, corresponding to the
+       boundary attributes marked in the array bdr_attr_is_ess.
+
+       For spaces with 'vdim' > 1, the 'component' array can be used
+       to restricts the marked tDOFs per boundary to the specified components.
+       If vdim > 1 then one can specify per boundary attribute which components
+       on a boundary are essential by assigning a value of true to its location
+       in the component array.
+       The component has dimensions number of boundary attributes x vdim. */
+   virtual void GetEssentialTrueDofs(const Array<int> &bdr_attr_is_ess,
+                                     Array<int> &ess_tdof_list,
+                                     const Array2D<bool> &component);
+
    /** If the given ldof is owned by the current processor, return its local
        tdof number, otherwise return -1 */
    int GetLocalTDofNumber(int ldof) const;

fem/pgridfunc.cpp

@@ -543,6 +543,12 @@ void ParGridFunction::ProjectCoefficient(Coefficient &coeff)
    }
 }
 
+//void ParGridFunction::ProjectBdrCoefficient(VectorFunctionCoefficient &vfcoeff,
+//                                            Array<int> &attr)
+//{
+//   GridFunction::ProjectBdrCoefficient(vfcoeff, attr);
+//}
+
 void ParGridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff)
 {
    // local maximal element attribute for each dof

fem/pgridfunc.hpp

@@ -238,6 +238,10 @@ public:
    using GridFunction::ProjectCoefficient;
    virtual void ProjectCoefficient(Coefficient &coeff);
 
+   //   using GridFunction::ProjectBdrCoefficient;
+   //   void ProjectBdrCoefficient(VectorFunctionCoefficient &vfcoeff,
+   //                              Array<int> &attr);
+
    using GridFunction::ProjectDiscCoefficient;
    /** @brief Project a discontinuous vector coefficient as a grid function on
        a continuous finite element space. The values in shared dofs are

fem/pnonlinearform.cpp

@@ -94,6 +94,195 @@ void ParNonlinearForm::Mult(const Vector &x, Vector &y) const
    MFEM_FORALL(i, N, Y_RW[idx[i]] = 0.0; );
 }
 
+void ParNonlinearForm::Mult(const Vector &x, Vector &y, ParMesh* pmesh,
+                            ParGridFunction* end_crds, ParGridFunction* beg_crds,
+                            const Vector &v) const
+{
+   //We're going to pretty much take everything from in NonlinearForm::Mult(x, y) into here
+   //with a few changes made to it
+
+   Array<int> vdofs;
+   Vector el_x, el_y, el_v;
+   const FiniteElement *fe;
+   IsoparametricTransformation beg_T;
+   IsoparametricTransformation end_T;
+   Mesh *mesh = fes->GetMesh();
+   const Vector &ptemp = Prolongate(x);
+   const Vector px(ptemp);
+   const Vector &ptemp1 = Prolongate(v);
+   const Vector pv(ptemp1);
+   Vector &py = P ? aux2.SetSize(P->Height()), aux2 : y;
+
+   py = 0.0;
+   GridFunction *nodes;
+
+   if (dnfi.Size())
+   {
+      for (int i = 0; i < fes->GetNE(); i++)
+      {
+         //Here we're pretty much forced to do things a little different
+         //from before since we have to explicitly call the pmesh GetElementTransformation
+         //if we don't then we end up with the exact same ElementTransformations
+         //for beg_T and end_T
+         fe = fes->GetFE(i);
+         fes->GetElementVDofs(i, vdofs);
+         pmesh->GetElementTransformation(i, &beg_T);
+         //These lines are updating
+         int own_nodes = 0;
+         nodes = end_crds;
+         pmesh->SwapNodes(nodes, own_nodes);
+         fe = fes->GetFE(i);
+         fes->GetElementVDofs(i, vdofs);
+         pmesh->GetElementTransformation(i, &end_T);
+         px.GetSubVector(vdofs, el_x);
+         pv.GetSubVector(vdofs, el_v);
+         for (int k = 0; k < dnfi.Size(); k++)
+         {
+            //This is the part that we're changing by adding more transformation terms
+            //as well
+            dnfi[k]->AssembleElementVector(*fe, beg_T, end_T, el_x, el_y, el_v);
+            py.AddElementVector(vdofs, el_y);
+         }
+         nodes = beg_crds;
+         pmesh->SwapNodes(nodes, own_nodes);
+      }
+   }
+   nodes = NULL;
+   //fix_me: Do we need to worry about the face transformations down here as well
+   //and what their element transformation is based upon? I'm feeling like the answer
+   //might be yes for certain applications. However, I'm not seeing the need for something
+   //like UMATs for the below at least currently.
+   if (fnfi.Size())
+   {
+      FaceElementTransformations *tr;
+      const FiniteElement *fe1, *fe2;
+      Array<int> vdofs2;
+
+      for (int i = 0; i < mesh->GetNumFaces(); i++)
+      {
+         tr = mesh->GetInteriorFaceTransformations(i);
+         if (tr != NULL)
+         {
+            fes->GetElementVDofs(tr->Elem1No, vdofs);
+            fes->GetElementVDofs(tr->Elem2No, vdofs2);
+            vdofs.Append (vdofs2);
+
+            px.GetSubVector(vdofs, el_x);
+
+            fe1 = fes->GetFE(tr->Elem1No);
+            fe2 = fes->GetFE(tr->Elem2No);
+
+            for (int k = 0; k < fnfi.Size(); k++)
+            {
+               fnfi[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
+               py.AddElementVector(vdofs, el_y);
+            }
+         }
+      }
+   }
+   //fix_me: Same as the previous fix_me in this function.
+   if (bfnfi.Size())
+   {
+      FaceElementTransformations *tr;
+      const FiniteElement *fe1, *fe2;
+
+      // Which boundary attributes need to be processed?
+      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++)
+      {
+         if (bfnfi_marker[k] == NULL)
+         {
+            bdr_attr_marker = 1;
+            break;
+         }
+         Array<int> &bdr_marker = *bfnfi_marker[k];
+         MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
+                     "invalid boundary marker for boundary face integrator #"
+                     << k << ", counting from zero");
+         for (int i = 0; i < bdr_attr_marker.Size(); i++)
+         {
+            bdr_attr_marker[i] |= bdr_marker[i];
+         }
+      }
+
+      for (int i = 0; i < fes -> GetNBE(); i++)
+      {
+         const int bdr_attr = mesh->GetBdrAttribute(i);
+         if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
+
+         tr = mesh->GetBdrFaceTransformations (i);
+         if (tr != NULL)
+         {
+            fes->GetElementVDofs(tr->Elem1No, vdofs);
+            px.GetSubVector(vdofs, el_x);
+
+            fe1 = fes->GetFE(tr->Elem1No);
+            // The fe2 object is really a dummy and not used on the boundaries,
+            // 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++)
+            {
+               if (bfnfi_marker[k] &&
+                   (*bfnfi_marker[k])[bdr_attr-1] == 0) { continue; }
+
+               bfnfi[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
+               py.AddElementVector(vdofs, el_y);
+            }
+         }
+      }
+   }
+
+   //The below is directly from ParNonlinearForm::Mult(x,y)
+
+   Y.SetData(aux2.GetData()); // aux2 contains A_local.P.x
+   //fix_me: same as the other previous fix_mes in this function
+   if (fnfi.Size())
+   {
+      // Terms over shared interior faces in parallel.
+      ParFiniteElementSpace *pfes = ParFESpace();
+      ParMesh *pmesh = pfes->GetParMesh();
+      FaceElementTransformations *tr;
+      const FiniteElement *fe1, *fe2;
+      Array<int> vdofs1, vdofs2;
+      Vector el_x, el_y;
+
+      X.SetData(aux1.GetData()); // aux1 contains P.x
+      X.ExchangeFaceNbrData();
+      const int n_shared_faces = pmesh->GetNSharedFaces();
+      for (int i = 0; i < n_shared_faces; i++)
+      {
+         tr = pmesh->GetSharedFaceTransformations(i, true);
+
+         fe1 = pfes->GetFE(tr->Elem1No);
+         fe2 = pfes->GetFaceNbrFE(tr->Elem2No);
+
+         pfes->GetElementVDofs(tr->Elem1No, vdofs1);
+         pfes->GetFaceNbrElementVDofs(tr->Elem2No, vdofs2);
+
+         el_x.SetSize(vdofs1.Size() + vdofs2.Size());
+         X.GetSubVector(vdofs1, el_x.GetData());
+         X.FaceNbrData().GetSubVector(vdofs2, el_x.GetData() + vdofs1.Size());
+
+         for (int k = 0; k < fnfi.Size(); k++)
+         {
+            fnfi[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
+            Y.AddElementVector(vdofs1, el_y.GetData());
+         }
+      }
+   }
+
+   P->MultTranspose(Y, y);
+
+   for (int i = 0; i < ess_tdof_list.Size(); i++)
+   {
+      y(ess_tdof_list[i]) = 0.0;
+   }
+}
+
+
 const SparseMatrix &ParNonlinearForm::GetLocalGradient(const Vector &x) const
 {
    MFEM_VERIFY(NonlinearForm::ext == nullptr,
@@ -104,6 +293,33 @@ const SparseMatrix &ParNonlinearForm::GetLocalGradient(const Vector &x) const
    return *Grad;
 }
 
+Operator &ParNonlinearForm::GetLocalGradient2(const Vector &x) const
+{
+   ParFiniteElementSpace *pfes = ParFESpace();
+
+   pGrad.Clear();
+
+   NonlinearForm::GetGradient(x); // (re)assemble Grad, no b.c.
+
+   OperatorHandle dA(pGrad.Type()), Ph(pGrad.Type());
+
+   if (fnfi.Size() == 0)
+   {
+      dA.MakeSquareBlockDiag(pfes->GetComm(), pfes->GlobalVSize(),
+                             pfes->GetDofOffsets(), Grad);
+   }
+   else
+   {
+      MFEM_ABORT("TODO: assemble contributions from shared face terms");
+   }
+
+   // TODO - construct Dof_TrueDof_Matrix directly in the pGrad format
+   Ph.ConvertFrom(pfes->Dof_TrueDof_Matrix());
+   pGrad.MakePtAP(dA, Ph);
+
+   return *pGrad.Ptr();
+}
+
 Operator &ParNonlinearForm::GetGradient(const Vector &x) const
 {
    if (NonlinearForm::ext) { return NonlinearForm::GetGradient(x); }

fem/pnonlinearform.hpp

@@ -18,6 +18,8 @@
 
 #include "pgridfunc.hpp"
 #include "nonlinearform.hpp"
+#include "../mesh/pmesh.hpp"
+
 
 namespace mfem
 {
@@ -51,10 +53,29 @@ public:
 
    virtual void Mult(const Vector &x, Vector &y) const;
 
+   /** The input essential dofs in @a x will, generally, be non-zero. However,
+       the output essential dofs in @a y will always be set to zero.
+       The input @a pmesh, @a end_crds @ beg_crds will be used to supply an element transformation
+       at a beggining time step and then one at an end time step. The ParMesh initially has its
+       nodes set at the beggining time step nodes and are updated to the current configuration location
+       before being returned to the beggining time step.
+       We're making sure the mesh returns to the same state we obtained it.
+       Both the input and the output vectors, @a x and @a y, must be true-dof
+       vectors, i.e. their size must be fes->GetTrueVSize().
+       The ParGridFunctions @a end_crds and @a beg_crds have virtual-dofs.*/
+   virtual void Mult(const Vector &x, Vector &y, ParMesh* pmesh,
+                     ParGridFunction* end_crds, ParGridFunction* beg_crds,
+                     const Vector &v) const;
+
+
    /// Return the local gradient matrix for the given true-dof vector x.
    /** The returned matrix does NOT have any boundary conditions imposed. */
    const SparseMatrix &GetLocalGradient(const Vector &x) const;
 
+   /// Return the local gradient matrix for the given true-dof vector x.
+   /** The returned matrix does NOT have any boundary conditions imposed. */
+   virtual Operator &GetLocalGradient2(const Vector &x) const;
+
    virtual Operator &GetGradient(const Vector &x) const;
 
    /// Set the operator type id for the parallel gradient matrix/operator.

linalg/vector.cpp

@@ -289,6 +289,30 @@ void Vector::SetVector(const Vector &v, int offset)
    }
 }
 
+void Vector::SetVector(const Vector &v, int offset, int sslen, int ssoffset)
+{
+   int vs = v.Size();
+   const double *vp = v.data + ssoffset;
+   double *p = data + offset;
+
+#ifdef MFEM_DEBUG
+   if (offset+sslen > size)
+   {
+      mfem_error("Vector::SetVector(const Vector &, int, int, int)");
+   }
+
+   if (ssoffset+sslen > vs)
+   {
+      mfem_error("Vector::SetVector(const Vector &, int, int, int)");
+   }
+#endif
+
+   for (int i = 0; i < sslen; i++)
+   {
+      p[i] = vp[i];
+   }
+}
+
 void Vector::Neg()
 {
    const bool use_dev = UseDevice();

linalg/vector.hpp

@@ -315,6 +315,11 @@ public:
 
    void SetVector(const Vector &v, int offset);
 
+   /// Set the subset of element in (*this) defined by @a offset
+   /// and the subset length, @a sslen, equal to the subset of
+   /// elements in @a v defined by @a ssoffset and @a sslen.
+   void SetVector (const Vector &v, int offset, int sslen, int ssoffset);
+
    /// (*this) = -(*this)
    void Neg();