MFEM  v4.0 Finite element discretization library
staticcond.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
10 // Software Foundation) version 2.1 dated February 1999.
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  double *A_data;
86  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:347
~StaticCondensation()
Destroy a StaticCondensation object.
Definition: staticcond.cpp:105
bool ReducesTrueVSize() const
Definition: staticcond.cpp:118
void ReduceRHS(const Vector &b, Vector &sc_b) const
Definition: staticcond.cpp:311
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:417
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:288
Abstract parallel finite element space.
Definition: pfespace.hpp:28
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:477
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:385
OpType * Is() const
Return the Operator pointer dynamically cast to a specified OpType.
Definition: handle.hpp:100
Data type sparse matrix.
Definition: sparsemat.hpp:40
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:350
Operator * Ptr() const
Access the underlying Operator pointer.
Definition: handle.hpp:82
Type
Enumeration defining IDs for some classes derived from Operator.
Definition: operator.hpp:135
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:85
void GetParallelMatrix(OperatorHandle &S_h) const
Return the parallel Schur complement matrix in the format specified by SetOperatorType().
Definition: staticcond.hpp:167
void Finalize()
Finalize the construction of the Schur complement matrix.
Definition: staticcond.cpp:235
void AssembleBdrMatrix(int el, const DenseMatrix &elmat)
Definition: staticcond.cpp:227
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:189
bool HasEliminatedBC() const
Return true if essential boundary conditions have been eliminated from the Schur complement matrix...
Definition: staticcond.hpp:142
static void ListToMarker(const Array< int > &list, int marker_size, Array< int > &marker, int mark_val=-1)
Definition: fespace.cpp:403
Vector data type.
Definition: vector.hpp:48
SparseMatrix & GetMatrix()
Return the serial Schur complement matrix.
Definition: staticcond.hpp:152
void ReduceSolution(const Vector &sol, Vector &sc_sol) const
Definition: staticcond.cpp:390
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:187
void Init(bool symmetric, bool block_diagonal)
Definition: staticcond.cpp:134
void ConvertListToReducedTrueDofs(const Array< int > &ess_tdof_list, Array< int > &ess_rtdof_list) const
Definition: staticcond.hpp:202
void ConvertMarkerToReducedTrueDofs(const Array< int > &ess_tdof_marker, Array< int > &ess_rtdof_marker) const
Definition: staticcond.cpp:440
void SetType(Operator::Type tid)
Invoke Clear() and set a new type id.
Definition: handle.hpp:124