MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
staticcond.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
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
22namespace 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 $$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}$$
48 where the block $A_1$ is itself block diagonal with small local blocks
49 and it is, therefore, easily invertible.
50
51 Starting with the block system
52 $$\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}$$
58 the reduced, statically condensed system is given by
59 $$S_{22} X_2 = B_2 - A_{21} A_{11}^{-1} B_1$$
60 where the Schur complement matrix $S_{22}$ is given by
61 $$S_{22} = A_{22} - A_{21} A_{11}^{-1} A_{12}.$$
62 After solving the Schur complement system, the $X_1$ part of the
63 solution can be recovered using the formula
64 $$X_1 = A_{11}^{-1} ( B_1 - A_{12} X_2 ).$$ */
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<real_t> A_data;
86 Memory<int> A_ipiv;
87
88 Array<int> ess_rtdof_list;
89
90public:
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.
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
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition fe_coll.hpp:27
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition fespace.hpp:716
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:667
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:648
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition fespace.hpp:713
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
Class used by MFEM to store pointers to host and/or device memory.
Pointer to an Operator of a specified type.
Definition handle.hpp:34
void SetType(Operator::Type tid)
Invoke Clear() and set a new type id.
Definition handle.hpp:132
Operator * Ptr() const
Access the underlying Operator pointer.
Definition handle.hpp:87
OpType * Is() const
Return the Operator pointer dynamically cast to a specified OpType.
Definition handle.hpp:108
DiagonalPolicy
Defines operator diagonal policy upon elimination of rows and/or columns.
Definition operator.hpp:48
Type
Enumeration defining IDs for some classes derived from Operator.
Definition operator.hpp:284
Abstract parallel finite element space.
Definition pfespace.hpp:29
Data type sparse matrix.
Definition sparsemat.hpp:51
HypreParMatrix & GetParallelMatrixElim()
Return the eliminated part of the parallel Schur complement matrix.
void AssembleMatrix(int el, const DenseMatrix &elmat)
int GetNPrDofs() const
Return the number of vector private dofs.
StaticCondensation(FiniteElementSpace *fespace)
Construct a StaticCondensation object.
void SetOperatorType(Operator::Type tid)
Set the operator type id for the parallel reduced matrix/operator.
void AssembleBdrMatrix(int el, const DenseMatrix &elmat)
void Finalize()
Finalize the construction of the Schur complement matrix.
void EliminateReducedTrueDofs(Matrix::DiagonalPolicy dpolicy)
Eliminate the internal reduced true dofs (set using SetEssentialTrueDofs()) from the Schur complement...
SparseMatrix & GetMatrix()
Return the serial Schur complement matrix.
void ConvertMarkerToReducedTrueDofs(const Array< int > &ess_tdof_marker, Array< int > &ess_rtdof_marker) const
void EliminateReducedTrueDofs(const Array< int > &ess_rtdof_list, Matrix::DiagonalPolicy dpolicy)
Eliminate the given reduced true dofs from the Schur complement matrix S.
bool ReducesTrueVSize() const
void Init(bool symmetric, bool block_diagonal)
void SetEssentialTrueDofs(const Array< int > &ess_tdof_list)
Determine and save internally essential reduced true dofs.
FiniteElementSpace * GetTraceFESpace()
Return a pointer to the reduced/trace FE space.
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....
void GetParallelMatrix(OperatorHandle &S_h) const
Return the parallel Schur complement matrix in the format specified by SetOperatorType().
void ComputeSolution(const Vector &b, const Vector &sc_sol, Vector &sol) const
HypreParMatrix & GetParallelMatrix()
Return the parallel Schur complement matrix.
void GetParallelMatrixElim(OperatorHandle &S_e_h) const
Return the eliminated part of the parallel Schur complement matrix in the format specified by SetOper...
ParFiniteElementSpace * GetParTraceFESpace()
Return a pointer to the parallel reduced/trace FE space.
SparseMatrix & GetMatrixElim()
Return the eliminated part of the serial Schur complement matrix.
void ReduceRHS(const Vector &b, Vector &sc_b) const
~StaticCondensation()
Destroy a StaticCondensation object.
void ConvertListToReducedTrueDofs(const Array< int > &ess_tdof_list_, Array< int > &ess_rtdof_list_) const
bool HasEliminatedBC() const
Return true if essential boundary conditions have been eliminated from the Schur complement matrix.
void ReduceSolution(const Vector &sol, Vector &sc_sol) const
int GetNExDofs() const
Return the number of vector exposed/reduced dofs.
Vector data type.
Definition vector.hpp:80
real_t b
Definition lissajous.cpp:42