Merge pull request mfem#4356 from mfem/najlkin/mixed-form-elim-dofs

Elimination of essential DOFs/BCs in MixedBilinearForm
adam-sim-dev · Sep 26, 2024 · 8ed11a5 · 8ed11a5
2 parents 0739640 + 0c91dba
commit 8ed11a5
Show file tree

Hide file tree

Showing 5 changed files with 132 additions and 53 deletions.
diff --git a/examples/ex36.cpp b/examples/ex36.cpp
@@ -258,7 +258,7 @@ int main(int argc, char *argv[])
          MixedBilinearForm a10(&H1fes,&L2fes);
          a10.AddDomainIntegrator(new MixedScalarMassIntegrator());
          a10.Assemble();
-         a10.EliminateTrialDofs(ess_bdr, x.GetBlock(0), rhs.GetBlock(1));
+         a10.EliminateTrialEssentialBC(ess_bdr, x.GetBlock(0), rhs.GetBlock(1));
          a10.Finalize();
          SparseMatrix &A10 = a10.SpMat();
 

diff --git a/examples/ex8.cpp b/examples/ex8.cpp
@@ -157,7 +157,7 @@ int main(int argc, char *argv[])
    MixedBilinearForm *B0 = new MixedBilinearForm(x0_space,test_space);
    B0->AddDomainIntegrator(new DiffusionIntegrator(one));
    B0->Assemble();
-   B0->EliminateTrialDofs(ess_bdr, x.GetBlock(x0_var), F);
+   B0->EliminateTrialEssentialBC(ess_bdr, x.GetBlock(x0_var), F);
    B0->Finalize();
 
    MixedBilinearForm *Bhat = new MixedBilinearForm(xhat_space,test_space);

diff --git a/fem/bilinearform.cpp b/fem/bilinearform.cpp
@@ -1941,36 +1941,59 @@ void MixedBilinearForm::AssembleBdrElementMatrix(
    mat->AddSubMatrix(test_vdofs_, trial_vdofs_, elmat, skip_zeros);
 }
 
-void MixedBilinearForm::EliminateTrialDofs (
+void MixedBilinearForm::EliminateTrialEssentialBC(
    const Array<int> &bdr_attr_is_ess, const Vector &sol, Vector &rhs )
 {
-   int i, j, k;
-   Array<int> tr_vdofs, cols_marker (trial_fes -> GetVSize());
+   Array<int> trial_ess_dofs;
+   trial_fes->GetEssentialVDofs(bdr_attr_is_ess, trial_ess_dofs);
+   mat->EliminateCols(trial_ess_dofs, &sol, &rhs);
+}
 
-   cols_marker = 0;
-   for (i = 0; i < trial_fes -> GetNBE(); i++)
-      if (bdr_attr_is_ess[trial_fes -> GetBdrAttribute (i)-1])
-      {
-         trial_fes -> GetBdrElementVDofs (i, tr_vdofs);
-         for (j = 0; j < tr_vdofs.Size(); j++)
-         {
-            if ( (k = tr_vdofs[j]) < 0 )
-            {
-               k = -1-k;
-            }
-            cols_marker[k] = 1;
-         }
-      }
-   mat -> EliminateCols (cols_marker, &sol, &rhs);
+void MixedBilinearForm::EliminateTrialEssentialBC(const Array<int>
+                                                  &bdr_attr_is_ess)
+{
+   Array<int> trial_ess_dofs;
+   trial_fes->GetEssentialVDofs(bdr_attr_is_ess, trial_ess_dofs);
+   mat->EliminateCols(trial_ess_dofs);
+}
+
+void MixedBilinearForm::EliminateTrialVDofs(const Array<int> &trial_vdofs_,
+                                            const Vector &sol, Vector &rhs)
+{
+   Array<int> trial_vdofs_marker;
+   FiniteElementSpace::ListToMarker(trial_vdofs_, mat->Width(),
+                                    trial_vdofs_marker);
+   mat->EliminateCols(trial_vdofs_marker, &sol, &rhs);
+}
+
+void MixedBilinearForm::EliminateTrialVDofs(const Array<int> &trial_vdofs_)
+{
+   if (mat_e == NULL)
+   {
+      mat_e = new SparseMatrix(mat->Height(), mat->Width());
+   }
+
+   Array<int> trial_vdofs_marker;
+   FiniteElementSpace::ListToMarker(trial_vdofs_, mat->Width(),
+                                    trial_vdofs_marker);
+   mat->EliminateCols(trial_vdofs_marker, *mat_e);
+   mat_e->Finalize();
 }
 
-void MixedBilinearForm::EliminateEssentialBCFromTrialDofs (
+void MixedBilinearForm::EliminateTrialVDofsInRHS(const Array<int> &trial_vdofs_,
+                                                 const Vector &x, Vector &b)
+{
+   mat_e->AddMult(x, b, -1.);
+}
+
+void MixedBilinearForm::EliminateEssentialBCFromTrialDofs(
    const Array<int> &marked_vdofs, const Vector &sol, Vector &rhs)
 {
-   mat -> EliminateCols (marked_vdofs, &sol, &rhs);
+   mat->EliminateCols(marked_vdofs, &sol, &rhs);
 }
 
-void MixedBilinearForm::EliminateTestDofs (const Array<int> &bdr_attr_is_ess)
+void MixedBilinearForm::EliminateTestEssentialBC(const Array<int>
+                                                 &bdr_attr_is_ess)
 {
    int i, j, k;
    Array<int> te_vdofs;
@@ -1990,6 +2013,14 @@ void MixedBilinearForm::EliminateTestDofs (const Array<int> &bdr_attr_is_ess)
       }
 }
 
+void MixedBilinearForm::EliminateTestVDofs(const Array<int> &test_vdofs_)
+{
+   for (int i=0; i<test_vdofs_.Size(); ++i)
+   {
+      mat->EliminateRow(test_vdofs_[i]);
+   }
+}
+
 void MixedBilinearForm::FormRectangularSystemMatrix(
    const Array<int> &trial_tdof_list,
    const Array<int> &test_tdof_list,
@@ -2026,20 +2057,9 @@ void MixedBilinearForm::FormRectangularSystemMatrix(
       mat = m;
    }
 
-   Array<int> ess_trial_tdof_marker, ess_test_tdof_marker;
-   FiniteElementSpace::ListToMarker(trial_tdof_list, trial_fes->GetTrueVSize(),
-                                    ess_trial_tdof_marker);
-   FiniteElementSpace::ListToMarker(test_tdof_list, test_fes->GetTrueVSize(),
-                                    ess_test_tdof_marker);
+   EliminateTrialVDofs(trial_tdof_list);
+   EliminateTestVDofs(test_tdof_list);
 
-   mat_e = new SparseMatrix(mat->Height(), mat->Width());
-   mat->EliminateCols(ess_trial_tdof_marker, *mat_e);
-
-   for (int i=0; i<test_tdof_list.Size(); ++i)
-   {
-      mat->EliminateRow(test_tdof_list[i]);
-   }
-   mat_e->Finalize();
    A.Reset(mat, false);
 }
 
@@ -2068,7 +2088,7 @@ void MixedBilinearForm::FormRectangularLinearSystem(
                                   A); // Set A = mat_e
    }
    // Eliminate essential BCs with B -= Ab xb
-   mat_e->AddMult(X, B, -1.0);
+   EliminateTrialVDofsInRHS(trial_tdof_list, X, B);
 
    B.SetSubVector(test_tdof_list, 0.0);
 }

diff --git a/fem/bilinearform.hpp b/fem/bilinearform.hpp
@@ -848,15 +848,37 @@ class MixedBilinearForm : public Matrix
    /** This will segfault if the usual sparse mat is not defined
        like when static condensation is being used or AllocMat() has
        not yet been called. */
-   const SparseMatrix &SpMat() const { return *mat; }
+   const SparseMatrix &SpMat() const
+   {
+      MFEM_VERIFY(mat, "mat is NULL and can't be dereferenced");
+      return *mat;
+   }
 
    /// Returns a reference to the sparse matrix:  $ M $
-   SparseMatrix &SpMat() { return *mat; }
+   SparseMatrix &SpMat()
+   {
+      MFEM_VERIFY(mat, "mat is NULL and can't be dereferenced");
+      return *mat;
+   }
 
    /**  @brief Nullifies the internal matrix $ M $ and returns a pointer
         to it.  Used for transferring ownership. */
    SparseMatrix *LoseMat() { SparseMatrix *tmp = mat; mat = NULL; return tmp; }
 
+   /// Returns a const reference to the sparse matrix of eliminated b.c.:  $ M_e $
+   const SparseMatrix &SpMatElim() const
+   {
+      MFEM_VERIFY(mat_e, "mat_e is NULL and can't be dereferenced");
+      return *mat_e;
+   }
+
+   /// Returns a reference to the sparse matrix of eliminated b.c.:  $ M_e $
+   SparseMatrix &SpMatElim()
+   {
+      MFEM_VERIFY(mat_e, "mat_e is NULL and can't be dereferenced");
+      return *mat_e;
+   }
+
    /// Adds a domain integrator. Assumes ownership of @a bfi.
    void AddDomainIntegrator(BilinearFormIntegrator *bfi);
 
@@ -1006,24 +1028,61 @@ class MixedBilinearForm : public Matrix
                                  Array<int> &test_vdofs,
                                  int skip_zeros = 1);
 
-   /// Eliminate essential boundary DOFs from the columns of the system.
+   /// Eliminate essential boundary trial DOFs from the system.
+   /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
+       the essential part of the boundary. */
+   void EliminateTrialEssentialBC(const Array<int> &bdr_attr_is_ess,
+                                  const Vector &sol, Vector &rhs);
+
+   /// Eliminate essential boundary trial DOFs from the system matrix.
    /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
-       the essential part of the boundary.  All entries in the columns will be
-       set to 0.0 through elimination.*/
-   void EliminateTrialDofs(const Array<int> &bdr_attr_is_ess,
-                           const Vector &sol, Vector &rhs);
-
-   /// Eliminate the list of DOFs from the columns of the system.
-   /** @a marked_vdofs is the of colunm numbers that will be eliminated.  All
-       entries in the columns will be set to 0.0 through elimination.*/
+       the essential part of the boundary. */
+   void EliminateTrialEssentialBC(const Array<int> &bdr_attr_is_ess);
+
+   /// (DEPRECATED) Eliminate essential boundary trial DOFs from the system.
+   /** @see EliminateTrialEssentialBC() */
+   MFEM_DEPRECATED void EliminateTrialDofs(const Array<int> &bdr_attr_is_ess,
+                                           const Vector &sol, Vector &rhs)
+   { EliminateTrialEssentialBC(bdr_attr_is_ess, sol, rhs); }
+
+   /// Eliminate the given trial @a vdofs. NOTE: here, @a vdofs is a list of DOFs.
+   /** In this case the eliminations are applied to the internal $ M $
+       and @a rhs without storing the elimination matrix $ M_e $. */
+   void EliminateTrialVDofs(const Array<int> &vdofs, const Vector &sol,
+                            Vector &rhs);
+
+   /// Eliminate the given trial @a vdofs, storing the eliminated part internally in $ M_e $.
+   /** This method works in conjunction with EliminateTrialVDofsInRHS() and allows
+       elimination of boundary conditions in multiple right-hand sides. In this
+       method, @a vdofs is a list of DOFs. */
+   void EliminateTrialVDofs(const Array<int> &vdofs);
+
+   /** @brief Use the stored eliminated part of the matrix (see
+       EliminateTrialVDofs(const Array<int> &)) to modify the r.h.s.
+       @a b; @a vdofs is a list of DOFs (non-directional, i.e. >= 0). */
+   void EliminateTrialVDofsInRHS(const Array<int> &vdofs, const Vector &x,
+                                 Vector &b);
+
+   /** @brief Similar to
+      EliminateTrialVDofs(const Array<int> &, const Vector &, Vector &)
+      but here @a ess_dofs is a marker (boolean) array on all vector-dofs
+      (@a ess_dofs[i] < 0 is true). */
    void EliminateEssentialBCFromTrialDofs(const Array<int> &marked_vdofs,
                                           const Vector &sol, Vector &rhs);
 
-   /// Eliminate essential boundary DOFs from the rows of the system.
+   /// Eliminate essential boundary test DOFs from the system matrix.
    /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
-       the essential part of the boundary.  All entries in the rows will be
-       set to 0.0 through elimination.*/
-   virtual void EliminateTestDofs(const Array<int> &bdr_attr_is_ess);
+       the essential part of the boundary. */
+   void EliminateTestEssentialBC(const Array<int> &bdr_attr_is_ess);
+
+   /// (DEPRECATED) Eliminate essential boundary test DOFs from the system.
+   /** @see EliminateTestEssentialBC() */
+   MFEM_DEPRECATED virtual void EliminateTestDofs(const Array<int>
+                                                  &bdr_attr_is_ess)
+   { EliminateTestEssentialBC(bdr_attr_is_ess); }
+
+   /// Eliminate the given test @a vdofs. NOTE: here, @a vdofs is a list of DOFs.
+   void EliminateTestVDofs(const Array<int> &vdofs);
 
    /** @brief Return in @a A that is column-constrained.
 

diff --git a/miniapps/solvers/block-solvers.cpp b/miniapps/solvers/block-solvers.cpp
@@ -168,7 +168,7 @@ DarcyProblem::DarcyProblem(Mesh &mesh, int num_refs, int order,
    bVarf_->AddDomainIntegrator(new VectorFEDivergenceIntegrator);
    bVarf_->Assemble();
    bVarf_->SpMat() *= -1.0;
-   bVarf_->EliminateTrialDofs(ess_bdr, u_, gform);
+   bVarf_->EliminateTrialEssentialBC(ess_bdr, u_, gform);
    bVarf_->Finalize();
    B_.Reset(bVarf_->ParallelAssemble());