MFEM  v4.4.0 Finite element discretization library
hybridization.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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_HYBRIDIZATION
13 #define MFEM_HYBRIDIZATION
14
15 #include "../config/config.hpp"
16 #include "fespace.hpp"
17 #include "bilininteg.hpp"
18
19 namespace mfem
20 {
21
22 /** Auxiliary class Hybridization, used to implement BilinearForm hybridization.
23
24  Hybridization can be viewed as a technique for solving linear systems
25  obtained through finite element assembly. The assembled matrix \f$A \f$ can
26  be written as:
27  \f[ A = P^T \hat{A} P, \f]
28  where \f$P \f$ is the matrix mapping the conforming finite element space to
29  the purely local finite element space without any inter-element constraints
30  imposed, and \f$\hat{A} \f$ is the block-diagonal matrix of all element
31  matrices.
32
33  We assume that:
34  - \f$\hat{A} \f$ is invertible,
35  - \f$P \f$ has a left inverse \f$R \f$, such that \f$R P = I \f$,
36  - a constraint matrix \f$C \f$ can be constructed, such that
37  \f$\operatorname{Ker}(C) = \operatorname{Im}(P) \f$.
38
39  Under these conditions, the linear system \f$A x = b \f$ can be solved
40  using the following procedure:
41  - solve for \f$\lambda \f$ in the linear system:
42  \f[ (C \hat{A}^{-1} C^T) \lambda = C \hat{A}^{-1} R^T b \f]
43  - compute \f$x = R \hat{A}^{-1} (R^T b - C^T \lambda) \f$
44
45  Hybridization is advantageous when the matrix
46  \f$H = (C \hat{A}^{-1} C^T) \f$ of the hybridized system is either smaller
47  than the original system, or is simpler to invert with a known method.
48
49  In some cases, e.g. high-order elements, the matrix \f$C \f$ can be written
50  as
51  \f[ C = \begin{pmatrix} 0 & C_b \end{pmatrix}, \f]
52  and then the hybridized matrix \f$H \f$ can be assembled using the identity
53  \f[ H = C_b S_b^{-1} C_b^T, \f]
54  where \f$S_b \f$ is the Schur complement of \f$\hat{A} \f$ with respect to
55  the same decomposition as the columns of \f$C \f$:
56  \f[ S_b = \hat{A}_b - \hat{A}_{bf} \hat{A}_{f}^{-1} \hat{A}_{fb}. \f]
57
58  Hybridization can also be viewed as a discretization method for imposing
59  (weak) continuity constraints between neighboring elements. */
61 {
62 protected:
65
67
70  double *Af_data;
71  int *Af_ipiv;
72
73 #ifdef MFEM_USE_MPI
74  HypreParMatrix *pC, *P_pc; // for parallel non-conforming meshes
76 #endif
77
78  void ConstructC();
79
80  void GetIBDofs(int el, Array<int> &i_dofs, Array<int> &b_dofs) const;
81
82  void GetBDofs(int el, int &num_idofs, Array<int> &b_dofs) const;
83
84  void ComputeH();
85
86  // Compute depending on mode:
87  // - mode 0: bf = Af^{-1} Rf^t b, where
88  // the non-"boundary" part of bf is set to 0;
89  // - mode 1: bf = Af^{-1} ( Rf^t b - Cf^t lambda ), where
90  // the "essential" part of bf is set to 0.
91  // Input: size(b) = fes->GetConformingVSize()
92  // size(lambda) = c_fes->GetConformingVSize()
93  void MultAfInv(const Vector &b, const Vector &lambda, Vector &bf,
94  int mode) const;
95
96 public:
97  /// Constructor
99  /// Destructor
100  ~Hybridization();
101
102  /** Set the integrator that will be used to construct the constraint matrix
103  C. The Hybridization object assumes ownership of the integrator, i.e. it
104  will delete the integrator when destroyed. */
106  { delete c_bfi; c_bfi = c_integ; }
107
108  /// Prepare the Hybridization object for assembly.
109  void Init(const Array<int> &ess_tdof_list);
110
111  /// Assemble the element matrix A into the hybridized system matrix.
112  void AssembleMatrix(int el, const DenseMatrix &A);
113
114  /// Assemble the boundary element matrix A into the hybridized system matrix.
115  void AssembleBdrMatrix(int bdr_el, const DenseMatrix &A);
116
117  /// Finalize the construction of the hybridized matrix.
118  void Finalize();
119
120  /// Return the serial hybridized matrix.
121  SparseMatrix &GetMatrix() { return *H; }
122
123 #ifdef MFEM_USE_MPI
124  /// Return the parallel hybridized matrix.
126
127  /** @brief Return the parallel hybridized matrix in the format specified by
128  SetOperatorType(). */
129  void GetParallelMatrix(OperatorHandle &H_h) const { H_h = pH; }
130
131  /// Set the operator type id for the parallel hybridized matrix/operator.
133 #endif
134
135  /** Perform the reduction of the given r.h.s. vector, b, to a r.h.s vector,
136  b_r, for the hybridized system. */
137  void ReduceRHS(const Vector &b, Vector &b_r) const;
138
139  /** Reconstruct the solution of the original system, sol, from solution of
140  the hybridized system, sol_r, and the original r.h.s. vector, b.
141  It is assumed that the vector sol has the right essential b.c. */
142  void ComputeSolution(const Vector &b, const Vector &sol_r,
143  Vector &sol) const;
144
145  /** @brief Destroy the current hybridization matrix while preserving the
146  computed constraint matrix and the set of essential true dofs. After
147  Reset(), a new hybridized matrix can be assembled via AssembleMatrix()
148  and Finalize(). The Mesh and FiniteElementSpace objects are assumed to be
149  un-modified. If that is not the case, a new Hybridization object must be
150  created. */
151  void Reset();
152 };
153
154 }
155
156 #endif
void GetParallelMatrix(OperatorHandle &H_h) const
Return the parallel hybridized matrix in the format specified by SetOperatorType().
Array< int > Af_offsets
FiniteElementSpace * c_fes
Array< int > Af_f_offsets
void SetConstraintIntegrator(BilinearFormIntegrator *c_integ)
void GetBDofs(int el, int &num_idofs, Array< int > &b_dofs) const
HypreParMatrix * P_pc
HypreParMatrix & GetParallelMatrix()
Return the parallel hybridized matrix.
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
void Reset()
Destroy the current hybridization matrix while preserving the computed constraint matrix and the set ...
SparseMatrix * Ct
void SetOperatorType(Operator::Type tid)
Set the operator type id for the parallel hybridized matrix/operator.
void Finalize()
Finalize the construction of the hybridized matrix.
OpType * Is() const
Return the Operator pointer dynamically cast to a specified OpType.
Definition: handle.hpp:108
void AssembleBdrMatrix(int bdr_el, const DenseMatrix &A)
Assemble the boundary element matrix A into the hybridized system matrix.
Data type sparse matrix.
Definition: sparsemat.hpp:46
SparseMatrix & GetMatrix()
Return the serial hybridized matrix.
double b
Definition: lissajous.cpp:42
void ComputeSolution(const Vector &b, const Vector &sol_r, Vector &sol) const
Array< int > hat_offsets
void MultAfInv(const Vector &b, const Vector &lambda, Vector &bf, int mode) const
BilinearFormIntegrator * c_bfi
Type
Enumeration defining IDs for some classes derived from Operator.
Definition: operator.hpp:255
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:34
SparseMatrix * H
void AssembleMatrix(int el, const DenseMatrix &A)
Assemble the element matrix A into the hybridized system matrix.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:88
void GetIBDofs(int el, Array< int > &i_dofs, Array< int > &b_dofs) const
OperatorHandle pH
~Hybridization()
Destructor.
FiniteElementSpace * fes
HypreParMatrix * pC
Array< int > hat_dofs_marker
Vector data type.
Definition: vector.hpp:60
void ReduceRHS(const Vector &b, Vector &b_r) const
void Init(const Array< int > &ess_tdof_list)
Prepare the Hybridization object for assembly.
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:337
Hybridization(FiniteElementSpace *fespace, FiniteElementSpace *c_fespace)
Constructor.
void SetType(Operator::Type tid)
Invoke Clear() and set a new type id.
Definition: handle.hpp:132