// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced // at the Lawrence Livermore National Laboratory. All Rights reserved. See files // LICENSE and NOTICE for details. LLNL-CODE-806117. // // This file is part of the MFEM library. For more information and source code // availability visit https://mfem.org. // // MFEM is free software; you can redistribute it and/or modify it under the // terms of the BSD-3 license. We welcome feedback and contributions, see file // CONTRIBUTING.md for details. #ifndef MFEM_HYBRIDIZATION #define MFEM_HYBRIDIZATION #include "../config/config.hpp" #include "fespace.hpp" #include "bilininteg.hpp" namespace mfem { /** Auxiliary class Hybridization, used to implement BilinearForm hybridization. Hybridization can be viewed as a technique for solving linear systems obtained through finite element assembly. The assembled matrix \f$ A \f$ can be written as: \f[ A = P^T \hat{A} P, \f] where \f$ P \f$ is the matrix mapping the conforming finite element space to the purely local finite element space without any inter-element constraints imposed, and \f$ \hat{A} \f$ is the block-diagonal matrix of all element matrices. We assume that: - \f$ \hat{A} \f$ is invertible, - \f$ P \f$ has a left inverse \f$ R \f$, such that \f$ R P = I \f$, - a constraint matrix \f$ C \f$ can be constructed, such that \f$ \operatorname{Ker}(C) = \operatorname{Im}(P) \f$. Under these conditions, the linear system \f$ A x = b \f$ can be solved using the following procedure: - solve for \f$ \lambda \f$ in the linear system: \f[ (C \hat{A}^{-1} C^T) \lambda = C \hat{A}^{-1} R^T b \f] - compute \f$ x = R \hat{A}^{-1} (R^T b - C^T \lambda) \f$ Hybridization is advantageous when the matrix \f$ H = (C \hat{A}^{-1} C^T) \f$ of the hybridized system is either smaller than the original system, or is simpler to invert with a known method. In some cases, e.g. high-order elements, the matrix \f$ C \f$ can be written as \f[ C = \begin{pmatrix} 0 & C_b \end{pmatrix}, \f] and then the hybridized matrix \f$ H \f$ can be assembled using the identity \f[ H = C_b S_b^{-1} C_b^T, \f] where \f$ S_b \f$ is the Schur complement of \f$ \hat{A} \f$ with respect to the same decomposition as the columns of \f$ C \f$: \f[ S_b = \hat{A}_b - \hat{A}_{bf} \hat{A}_{f}^{-1} \hat{A}_{fb}. \f] Hybridization can also be viewed as a discretization method for imposing (weak) continuity constraints between neighboring elements. */ class Hybridization { protected: FiniteElementSpace *fes, *c_fes; BilinearFormIntegrator *c_bfi; SparseMatrix *Ct, *H; Array<int> hat_offsets, hat_dofs_marker; Array<int> Af_offsets, Af_f_offsets; double *Af_data; int *Af_ipiv; #ifdef MFEM_USE_MPI HypreParMatrix *pC, *P_pc; // for parallel non-conforming meshes OperatorHandle pH; #endif void ConstructC(); void GetIBDofs(int el, Array<int> &i_dofs, Array<int> &b_dofs) const; void GetBDofs(int el, int &num_idofs, Array<int> &b_dofs) const; void ComputeH(); // Compute depending on mode: // - mode 0: bf = Af^{-1} Rf^t b, where // the non-"boundary" part of bf is set to 0; // - mode 1: bf = Af^{-1} ( Rf^t b - Cf^t lambda ), where // the "essential" part of bf is set to 0. // Input: size(b) = fes->GetConformingVSize() // size(lambda) = c_fes->GetConformingVSize() void MultAfInv(const Vector &b, const Vector &lambda, Vector &bf, int mode) const; public: /// Constructor Hybridization(FiniteElementSpace *fespace, FiniteElementSpace *c_fespace); /// Destructor ~Hybridization(); /** Set the integrator that will be used to construct the constraint matrix C. The Hybridization object assumes ownership of the integrator, i.e. it will delete the integrator when destroyed. */ void SetConstraintIntegrator(BilinearFormIntegrator *c_integ) { delete c_bfi; c_bfi = c_integ; } /// Prepare the Hybridization object for assembly. void Init(const Array<int> &ess_tdof_list); /// Assemble the element matrix A into the hybridized system matrix. void AssembleMatrix(int el, const DenseMatrix &A); /// Assemble the boundary element matrix A into the hybridized system matrix. void AssembleBdrMatrix(int bdr_el, const DenseMatrix &A); /// Finalize the construction of the hybridized matrix. void Finalize(); /// Return the serial hybridized matrix. SparseMatrix &GetMatrix() { return *H; } #ifdef MFEM_USE_MPI /// Return the parallel hybridized matrix. HypreParMatrix &GetParallelMatrix() { return *pH.Is<HypreParMatrix>(); } /** @brief Return the parallel hybridized matrix in the format specified by SetOperatorType(). */ void GetParallelMatrix(OperatorHandle &H_h) const { H_h = pH; } /// Set the operator type id for the parallel hybridized matrix/operator. void SetOperatorType(Operator::Type tid) { pH.SetType(tid); } #endif /** Perform the reduction of the given r.h.s. vector, b, to a r.h.s vector, b_r, for the hybridized system. */ void ReduceRHS(const Vector &b, Vector &b_r) const; /** Reconstruct the solution of the original system, sol, from solution of the hybridized system, sol_r, and the original r.h.s. vector, b. It is assumed that the vector sol has the right essential b.c. */ void ComputeSolution(const Vector &b, const Vector &sol_r, Vector &sol) const; /** @brief Destroy the current hybridization matrix while preserving the computed constraint matrix and the set of essential true dofs. After Reset(), a new hybridized matrix can be assembled via AssembleMatrix() and Finalize(). The Mesh and FiniteElementSpace objects are assumed to be un-modified. If that is not the case, a new Hybridization object must be created. */ void Reset(); }; } #endif