MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
staticcond.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #ifndef MFEM_STATIC_CONDENSATION
13 #define MFEM_STATIC_CONDENSATION
14 
15 #include "../config/config.hpp"
16 #include "fespace.hpp"
17 
18 #ifdef MFEM_USE_MPI
19 #include "pfespace.hpp"
20 #endif
21 
22 namespace mfem
23 {
24 
25 /** Auxiliary class StaticCondensation, used to implement static condensation
26  in class BilinearForm.
27 
28  Static condensation is a technique for solving linear systems by eliminating
29  groups/blocks of unknowns and reducing the original system to the remaining
30  interfacial unknowns. The assumption is that unknowns in one group are
31  connected (in the graph of the matrix) only to unknowns in the same group
32  or to interfacial unknowns but not to other groups.
33 
34  For finite element systems, the groups correspond to degrees of freedom
35  (DOFs) associated with the interior of the elements. The rest of the DOFs
36  (associated with the element boundaries) are interfacial.
37 
38  In block form the matrix of the system can be written as
39  \f[ A =
40  \begin{pmatrix}
41  A_{11} & A_{12} \\
42  A_{21} & A_{22}
43  \end{pmatrix}
44  \begin{array}{l}
45  \text{-- groups: element interior/private DOFs} \\
46  \text{-- interface: element boundary/exposed DOFs}
47  \end{array} \f]
48  where the block \f$ A_1 \f$ is itself block diagonal with small local blocks
49  and it is, therefore, easily invertible.
50 
51  Starting with the block system
52  \f[ \begin{pmatrix}
53  A_{11} & A_{12} \\
54  A_{21} & A_{22}
55  \end{pmatrix}
56  \begin{pmatrix} X_1 \\ X_2 \end{pmatrix} =
57  \begin{pmatrix} B_1 \\ B_2 \end{pmatrix} \f]
58  the reduced, statically condensed system is given by
59  \f[ S_{22} X_2 = B_2 - A_{21} A_{11}^{-1} B_1 \f]
60  where the Schur complement matrix \f$ S_{22} \f$ is given by
61  \f[ S_{22} = A_{22} - A_{21} A_{11}^{-1} A_{12}. \f]
62  After solving the Schur complement system, the \f$ X_1 \f$ part of the
63  solution can be recovered using the formula
64  \f[ X_1 = A_{11}^{-1} ( B_1 - A_{12} X_2 ). \f] */
66 {
67  FiniteElementSpace *fes, *tr_fes;
69  Table elem_pdof; // Element to private dof
70  int npdofs; // Number of private dofs
71  Array<int> rdof_edof; // Map from reduced dofs to exposed dofs
72 
73  // Schur complement: S = A_ee - A_ep (A_pp)^{-1} A_pe.
74  SparseMatrix *S, *S_e;
75 #ifdef MFEM_USE_MPI
76  ParFiniteElementSpace *pfes, *tr_pfes;
77  OperatorHandle pS, pS_e;
78  bool Parallel() const { return (tr_pfes != NULL); }
79 #else
80  bool Parallel() const { return false; }
81 #endif
82 
83  bool symm; // TODO: handle the symmetric case correctly.
84  Array<int> A_offsets, A_ipiv_offsets;
85  Memory<double> A_data;
86  Memory<int> A_ipiv;
87 
88  Array<int> ess_rtdof_list;
89 
90 public:
91  /// Construct a StaticCondensation object.
93  /// Destroy a StaticCondensation object.
95 
96  /// Return the number of vector private dofs.
97  int GetNPrDofs() const { return npdofs; }
98  /// Return the number of vector exposed/reduced dofs.
99  int GetNExDofs() const { return tr_fes->GetVSize(); }
100  /** Return true if applying the static condensation actually reduces the
101  (global) number of true vector dofs. */
102  bool ReducesTrueVSize() const;
103 
104  /** Prepare the StaticCondensation object to assembly: allocate the Schur
105  complement matrix and the other element-wise blocks. */
106  void Init(bool symmetric, bool block_diagonal);
107 
108  /// Return a pointer to the reduced/trace FE space.
109  FiniteElementSpace *GetTraceFESpace() { return tr_fes; }
110 
111 #ifdef MFEM_USE_MPI
112  /// Return a pointer to the parallel reduced/trace FE space.
114 #endif
115  /** Assemble the contribution to the Schur complement from the given
116  element matrix 'elmat'; save the other blocks internally: A_pp_inv, A_pe,
117  and A_ep. */
118  void AssembleMatrix(int el, const DenseMatrix &elmat);
119 
120  /** Assemble the contribution to the Schur complement from the given boundary
121  element matrix 'elmat'. */
122  void AssembleBdrMatrix(int el, const DenseMatrix &elmat);
123 
124  /// Finalize the construction of the Schur complement matrix.
125  void Finalize();
126 
127  /// Determine and save internally essential reduced true dofs.
128  void SetEssentialTrueDofs(const Array<int> &ess_tdof_list)
129  { ConvertListToReducedTrueDofs(ess_tdof_list, ess_rtdof_list); }
130 
131  /// Eliminate the given reduced true dofs from the Schur complement matrix S.
132  void EliminateReducedTrueDofs(const Array<int> &ess_rtdof_list,
133  Matrix::DiagonalPolicy dpolicy);
134 
135  /// @brief Eliminate the internal reduced true dofs (set using
136  /// SetEssentialTrueDofs()) from the Schur complement matrix S.
138  { EliminateReducedTrueDofs(ess_rtdof_list, dpolicy); }
139 
140  /** @brief Return true if essential boundary conditions have been eliminated
141  from the Schur complement matrix. */
142  bool HasEliminatedBC() const
143  {
144 #ifndef MFEM_USE_MPI
145  return S_e;
146 #else
147  return S_e || pS_e.Ptr();
148 #endif
149  }
150 
151  /// Return the serial Schur complement matrix.
152  SparseMatrix &GetMatrix() { return *S; }
153 
154  /// Return the eliminated part of the serial Schur complement matrix.
155  SparseMatrix &GetMatrixElim() { return *S_e; }
156 
157 #ifdef MFEM_USE_MPI
158  /// Return the parallel Schur complement matrix.
160 
161  /// Return the eliminated part of the parallel Schur complement matrix.
163  { return *pS_e.Is<HypreParMatrix>(); }
164 
165  /** @brief Return the parallel Schur complement matrix in the format
166  specified by SetOperatorType(). */
167  void GetParallelMatrix(OperatorHandle &S_h) const { S_h = pS; }
168 
169  /** @brief Return the eliminated part of the parallel Schur complement matrix
170  in the format specified by SetOperatorType(). */
171  void GetParallelMatrixElim(OperatorHandle &S_e_h) const { S_e_h = pS_e; }
172 
173  /// Set the operator type id for the parallel reduced matrix/operator.
175  { pS.SetType(tid); pS_e.SetType(tid); }
176 #endif
177 
178  /** Given a RHS vector for the full linear system, compute the RHS for the
179  reduced linear system: sc_b = b_e - A_ep A_pp_inv b_p. */
180  void ReduceRHS(const Vector &b, Vector &sc_b) const;
181 
182  /** Restrict a solution vector on the full FE space dofs to a vector on the
183  reduced/trace true FE space dofs. */
184  void ReduceSolution(const Vector &sol, Vector &sc_sol) const;
185 
186  /** @brief Set the reduced solution `X` and r.h.s `B` vectors from the full
187  linear system solution `x` and r.h.s. `b` vectors.
188 
189  This method should be called after the internal reduced essential dofs
190  have been set using SetEssentialTrueDofs() and both the Schur complement
191  and its eliminated part have been finalized. */
192  void ReduceSystem(Vector &x, Vector &b, Vector &X, Vector &B,
193  int copy_interior = 0) const;
194 
195  /** Restrict a marker Array on the true FE space dofs to a marker Array on
196  the reduced/trace true FE space dofs. */
197  void ConvertMarkerToReducedTrueDofs(const Array<int> &ess_tdof_marker,
198  Array<int> &ess_rtdof_marker) const;
199 
200  /** Restrict a list of true FE space dofs to a list of reduced/trace true FE
201  space dofs. */
202  void ConvertListToReducedTrueDofs(const Array<int> &ess_tdof_list_,
203  Array<int> &ess_rtdof_list_) const
204  {
205  Array<int> ess_tdof_marker, ess_rtdof_marker;
206  FiniteElementSpace::ListToMarker(ess_tdof_list_, fes->GetTrueVSize(),
207  ess_tdof_marker);
208  ConvertMarkerToReducedTrueDofs(ess_tdof_marker, ess_rtdof_marker);
209  FiniteElementSpace::MarkerToList(ess_rtdof_marker, ess_rtdof_list_);
210  }
211 
212  /** Given a solution of the reduced system 'sc_sol' and the RHS 'b' for the
213  full linear system, compute the solution of the full system 'sol'. */
214  void ComputeSolution(const Vector &b, const Vector &sc_sol,
215  Vector &sol) const;
216 };
217 
218 }
219 
220 #endif
void SetOperatorType(Operator::Type tid)
Set the operator type id for the parallel reduced matrix/operator.
Definition: staticcond.hpp:174
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:563
~StaticCondensation()
Destroy a StaticCondensation object.
Definition: staticcond.cpp:103
bool ReducesTrueVSize() const
Definition: staticcond.cpp:116
void ReduceRHS(const Vector &b, Vector &sc_b) const
Definition: staticcond.cpp:309
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
ParFiniteElementSpace * GetParTraceFESpace()
Return a pointer to the parallel reduced/trace FE space.
Definition: staticcond.hpp:113
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
void ReduceSystem(Vector &x, Vector &b, Vector &X, Vector &B, int copy_interior=0) const
Set the reduced solution X and r.h.s B vectors from the full linear system solution x and r...
Definition: staticcond.cpp:416
void EliminateReducedTrueDofs(const Array< int > &ess_rtdof_list, Matrix::DiagonalPolicy dpolicy)
Eliminate the given reduced true dofs from the Schur complement matrix S.
Definition: staticcond.cpp:286
Abstract parallel finite element space.
Definition: pfespace.hpp:28
void ConvertListToReducedTrueDofs(const Array< int > &ess_tdof_list_, Array< int > &ess_rtdof_list_) const
Definition: staticcond.hpp:202
void GetParallelMatrixElim(OperatorHandle &S_e_h) const
Return the eliminated part of the parallel Schur complement matrix in the format specified by SetOper...
Definition: staticcond.hpp:171
StaticCondensation(FiniteElementSpace *fespace)
Construct a StaticCondensation object.
Definition: staticcond.cpp:17
void ComputeSolution(const Vector &b, const Vector &sc_sol, Vector &sol) const
Definition: staticcond.cpp:476
int GetNExDofs() const
Return the number of vector exposed/reduced dofs.
Definition: staticcond.hpp:99
static void MarkerToList(const Array< int > &marker, Array< int > &list)
Convert a Boolean marker array to a list containing all marked indices.
Definition: fespace.cpp:598
OpType * Is() const
Return the Operator pointer dynamically cast to a specified OpType.
Definition: handle.hpp:108
Data type sparse matrix.
Definition: sparsemat.hpp:46
double b
Definition: lissajous.cpp:42
HypreParMatrix & GetParallelMatrix()
Return the parallel Schur complement matrix.
Definition: staticcond.hpp:159
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:566
Operator * Ptr() const
Access the underlying Operator pointer.
Definition: handle.hpp:87
Type
Enumeration defining IDs for some classes derived from Operator.
Definition: operator.hpp:255
int GetNPrDofs() const
Return the number of vector private dofs.
Definition: staticcond.hpp:97
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:88
void GetParallelMatrix(OperatorHandle &S_h) const
Return the parallel Schur complement matrix in the format specified by SetOperatorType().
Definition: staticcond.hpp:167
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
void Finalize()
Finalize the construction of the Schur complement matrix.
Definition: staticcond.cpp:233
void AssembleBdrMatrix(int el, const DenseMatrix &elmat)
Definition: staticcond.cpp:225
HypreParMatrix & GetParallelMatrixElim()
Return the eliminated part of the parallel Schur complement matrix.
Definition: staticcond.hpp:162
FiniteElementSpace * GetTraceFESpace()
Return a pointer to the reduced/trace FE space.
Definition: staticcond.hpp:109
void AssembleMatrix(int el, const DenseMatrix &elmat)
Definition: staticcond.cpp:187
bool HasEliminatedBC() const
Return true if essential boundary conditions have been eliminated from the Schur complement matrix...
Definition: staticcond.hpp:142
DiagonalPolicy
Defines operator diagonal policy upon elimination of rows and/or columns.
Definition: operator.hpp:46
static void ListToMarker(const Array< int > &list, int marker_size, Array< int > &marker, int mark_val=-1)
Convert an array of indices (list) to a Boolean marker array where all indices in the list are marked...
Definition: fespace.cpp:617
Vector data type.
Definition: vector.hpp:60
SparseMatrix & GetMatrix()
Return the serial Schur complement matrix.
Definition: staticcond.hpp:152
void ReduceSolution(const Vector &sol, Vector &sc_sol) const
Definition: staticcond.cpp:389
SparseMatrix & GetMatrixElim()
Return the eliminated part of the serial Schur complement matrix.
Definition: staticcond.hpp:155
void SetEssentialTrueDofs(const Array< int > &ess_tdof_list)
Determine and save internally essential reduced true dofs.
Definition: staticcond.hpp:128
void EliminateReducedTrueDofs(Matrix::DiagonalPolicy dpolicy)
Eliminate the internal reduced true dofs (set using SetEssentialTrueDofs()) from the Schur complement...
Definition: staticcond.hpp:137
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:337
void Init(bool symmetric, bool block_diagonal)
Definition: staticcond.cpp:132
void ConvertMarkerToReducedTrueDofs(const Array< int > &ess_tdof_marker, Array< int > &ess_rtdof_marker) const
Definition: staticcond.cpp:439
void SetType(Operator::Type tid)
Invoke Clear() and set a new type id.
Definition: handle.hpp:132