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