Merge branch 'master' into feature/partial_ess_bcs

Claude helped generate the docs / do the refactor as I was lazy...
However, the unit tests pass and the amount of common code between the two cases has now been reduced which is always a good thing...

fem/fespace.cpp

@@ -547,11 +547,22 @@ void MarkDofs(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> 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++)
    {
       if (bdr_attr_is_ess[GetBdrAttribute(i)-1])
@@ -663,6 +674,54 @@ void FiniteElementSpace::GetEssentialTrueDofs(const Array<int> &bdr_attr_is_ess,
    MarkerToList(ess_tdofs, ess_tdof_list);
 }
 
+void FiniteElementSpace::GetEssentialVDofsFromComponent(
+   const Array<int> &bdr_attr_is_ess,
+   const Array2D<bool> &component,
+   Array<int> &ess_vdofs) const
+{
+   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> 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++)
+   {
+      const bool overwrite = (i == 0);
+      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);
+   }
+}
+
+void FiniteElementSpace::GetEssentialTrueDofs(const Array<int> &bdr_attr_is_ess,
+                                              Array<int> &ess_tdof_list,
+                                              const Array2D<bool> &component)
+{
+   Array<int> ess_vdofs, ess_tdofs;
+
+   GetEssentialVDofsFromComponent(bdr_attr_is_ess, component, ess_vdofs);
+
+   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

@@ -595,6 +595,29 @@ protected:
    virtual void CopyProlongationAndRestriction(const FiniteElementSpace &fes,
                                                const Array<int> *perm);
 
+   /** @brief Helper function to mark essential VDofs based on a component
+       matrix specifying which vector components are essential on each boundary
+       attribute.
+
+       This method loops over all vector dimensions (vdim) and calls
+       GetEssentialVDofs() for each component where the corresponding entry
+       in the @a component matrix is true.
+
+       @param[in] bdr_attr_is_ess Array marking which boundary attributes are
+                                  essential (1 for essential, 0 otherwise).
+       @param[in] component 2D boolean array of size (num_bdr_attributes x vdim)
+                           indicating which vector components are essential for
+                           each boundary attribute. component(j,i) == true means
+                           component i is essential on boundary attribute j.
+       @param[out] ess_vdofs Marker array for essential VDofs. On exit,
+                            ess_vdofs[i] != 0 if VDof i is essential.
+
+       @note This is a helper method used by GetEssentialTrueDofs() when
+             component-wise boundary conditions are specified. */
+   void GetEssentialVDofsFromComponent(const Array<int> &bdr_attr_is_ess,
+                                       const Array2D<bool> &component,
+                                       Array<int> &ess_vdofs) const;
+
 public:
 
 
@@ -1356,11 +1379,17 @@ 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 which allows the accumulation of multiple DOFs into a single array.
+       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.
@@ -1370,6 +1399,20 @@ public:
                                      Array<int> &ess_tdof_list,
                                      int component = -1) 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.
+
+       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/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

@@ -165,14 +165,39 @@ public:
    const Array<NonlinearFormIntegrator*> &GetBdrFaceIntegrators() const
    { return bfnfi; }
 
-   /// Specify essential boundary conditions.
-   /** This method calls FiniteElementSpace::GetEssentialTrueDofs() and stores
+   /** @brief Specify essential boundary conditions.
+
+       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. */
+       have zero entries at the essential true dofs.
+
+       @note The values in the essential vdofs have to come from the initial guess.
+       */
    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.
+
+       @note The values in the essential vdofs have to come from the initial guess.
+       */
+   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. */
    void SetEssentialVDofs(const Array<int> &ess_vdofs_list);

fem/pfespace.cpp

@@ -1122,9 +1122,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).
@@ -1184,6 +1186,39 @@ 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_VERIFY(!IsVariableOrderH1(),
+               "Variable order H1 spaces are currently not supported with this feature");
+
+   Array<int> ess_dofs, true_ess_dofs;
+
+   GetEssentialVDofsFromComponent(bdr_attr_is_ess, component, ess_dofs);
+   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);
+}
+
 void ParFiniteElementSpace::GetEssentialTrueDofsVar(const Array<int>
                                                     &bdr_attr_is_ess,
                                                     const Array<int> &ess_dofs,

fem/pfespace.hpp

@@ -420,10 +420,18 @@ public:
        "partially conforming") space. */
    void Synchronize(Array<int> &ldof_marker) const;
 
-   /// Determine the boundary degrees of freedom
-   void GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
-                          Array<int> &ess_dofs,
-                          int component = -1) const override;
+   /** @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,
+                                  bool overwrite = true) const override;
 
    /** Get a list of essential true dofs, ess_tdof_list, corresponding to the
        boundary attributes marked in the array bdr_attr_is_ess. */
@@ -448,6 +456,19 @@ public:
    void GetExteriorTrueDofs(Array<int> &ext_tdof_list,
                             int component = -1) const override;
 
+   /** @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) override;
+
    /** If the given ldof is owned by the current processor, return its local
        tdof number, otherwise return -1 */
    int GetLocalTDofNumber(int ldof) const;

tests/unit/CMakeLists.txt

@@ -122,6 +122,7 @@ set(UNIT_TESTS_SRCS
   fem/test_fe_pos.cpp
   fem/test_fe_symmetry.cpp
   fem/test_fe.cpp
+  fem/test_fespace_get_ess_true_dofs.cpp
   fem/test_get_value.cpp
   fem/test_getderivative.cpp
   fem/test_getgradient.cpp

tests/unit/fem/test_fespace_get_ess_true_dofs.cpp

@@ -0,0 +1,74 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+using namespace mfem;
+
+TEST_CASE("FESpace Get Essential True DOFs",
+          "[FESpace Get Essential True DOFs]")
+{
+   std::cout << "Testing get essential true dofs" << std::endl;
+   int order_h1 = 3, n = 2, dim = 3;
+
+   Mesh mesh = Mesh::MakeCartesian3D(
+                  n, n, n, Element::HEXAHEDRON, 1.0, 1.0, 1.0);
+   mesh.SetCurvature(order_h1);
+
+   H1_FECollection fec(order_h1, dim);
+   FiniteElementSpace fe_space(&mesh, &fec, dim);
+
+
+   const int num_bdr_attr = fe_space.GetMesh()->bdr_attributes.Max();
+
+
+   Array<int> ess_tdofs_2d, ess_tdofs_1d, ess_tdofs_tmp, ess_bdrs;
+   Array2D<bool> comps(num_bdr_attr, dim);
+   ess_bdrs.SetSize(num_bdr_attr);
+   comps = false;
+   ess_bdrs = 0;
+
+   // simple xy boundary condition on all surfaces
+   // could do something more complex but don't really want to...
+   for (int i = 0; i < num_bdr_attr; i++)
+   {
+      ess_bdrs[i] = 1;
+      comps(i, 0) = true;
+      comps(i, 2) = true;
+   }
+
+   fe_space.GetEssentialTrueDofs(ess_bdrs, ess_tdofs_2d, comps);
+
+   // Now for the old way
+   fe_space.GetEssentialTrueDofs(ess_bdrs, ess_tdofs_tmp, 0);
+   ess_tdofs_1d.Append(ess_tdofs_tmp);
+   ess_tdofs_tmp.DeleteAll();
+   fe_space.GetEssentialTrueDofs(ess_bdrs, ess_tdofs_tmp, 2);
+   ess_tdofs_1d.Append(ess_tdofs_tmp);
+   // Sort the 2 arrays in order to compare them
+   ess_tdofs_2d.Sort();
+   ess_tdofs_1d.Sort();
+
+   int diff = 0;
+
+   for (int i = 0; i < ess_tdofs_2d.Size(); i++)
+   {
+      diff += std::abs(ess_tdofs_2d[i] - ess_tdofs_1d[i]);
+   }
+
+   std::cout << "Difference in essential tdofs approaches is: " << diff <<
+             std::endl;
+
+   REQUIRE(diff == 0);
+
+}
+