1 // Copyright (c) 2010-2021, 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_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. */ 60 class Hybridization 61 { 62 protected: 63 FiniteElementSpace *fes, *c_fes; 64 BilinearFormIntegrator *c_bfi; 65 66 SparseMatrix *Ct, *H; 67 68 Array<int> hat_offsets, hat_dofs_marker; 69 Array<int> Af_offsets, Af_f_offsets; 70 double *Af_data; 71 int *Af_ipiv; 72 73 #ifdef MFEM_USE_MPI 74 HypreParMatrix *pC, *P_pc; // for parallel non-conforming meshes 75 OperatorHandle pH; 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 98 Hybridization(FiniteElementSpace *fespace, FiniteElementSpace *c_fespace); 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. */ SetConstraintIntegrator(BilinearFormIntegrator * c_integ)105 void SetConstraintIntegrator(BilinearFormIntegrator *c_integ) 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. GetMatrix()121 SparseMatrix &GetMatrix() { return *H; } 122 123 #ifdef MFEM_USE_MPI 124 /// Return the parallel hybridized matrix. GetParallelMatrix()125 HypreParMatrix &GetParallelMatrix() { return *pH.Is<HypreParMatrix>(); } 126 127 /** @brief Return the parallel hybridized matrix in the format specified by 128 SetOperatorType(). */ GetParallelMatrix(OperatorHandle & H_h) const129 void GetParallelMatrix(OperatorHandle &H_h) const { H_h = pH; } 130 131 /// Set the operator type id for the parallel hybridized matrix/operator. SetOperatorType(Operator::Type tid)132 void SetOperatorType(Operator::Type tid) { pH.SetType(tid); } 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 157