// 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_BLOCK_STATIC_CONDENSATION #define MFEM_BLOCK_STATIC_CONDENSATION #include "mfem.hpp" namespace mfem { /** @brief Class that performs static condensation of interior dofs for multiple FE spaces. This class is used in class DPGWeakFrom. It is suitable for systems resulting from the discretization of multiple FE spaces. It eliminates the dofs associated with the interior of the elements and returns the reduced system which contains only the interfacial dofs. The ordering of the dofs in the matrix is implied by the ordering given by the FE spaces, but there is no assumption on the ordering of the FE spaces. This class handles both serial and parallel FE spaces. */ class BlockStaticCondensation { int height, width; int nblocks; // original number of blocks int rblocks; // reduces number of blocks Mesh * mesh = nullptr; bool parallel = false; // original set of Finite Element Spaces Array fes; // indicates if the original space is a trace space Array IsTraceSpace; // New set of "reduced" Finite Element Spaces // (after static condensation) Array tr_fes; Array dof_offsets; Array tdof_offsets; Array rdof_offsets; Array rtdof_offsets; /** Schur complement matrix S = A_ii - A_ib (A_bb)^{-1} A_bi. Here (b) and (i) stand for "bubble" and "interface" respectively */ BlockMatrix * S = nullptr; BlockMatrix * S_e = nullptr; BlockVector * y = nullptr; Array lmat; Array lvec; Array rdof_edof; // Map from reduced dofs to exposed dofs Array ess_rtdof_list; BlockMatrix * P = nullptr; // Block Prolongation BlockMatrix * R = nullptr; // Block Restriction #ifdef MFEM_USE_MPI BlockOperator * pS = nullptr; BlockOperator * pS_e = nullptr; BlockOperator * pP = nullptr; #endif bool Parallel() const { return parallel; } // tr_idx (trace dofs indices) // int_idx (interior dof indices) void GetReducedElementIndicesAndOffsets(int el, Array & tr_idx, Array & int_idx, Array & offsets) const; void GetReducedElementVDofs(int el, Array & rdofs) const; void GetElementVDofs(int el, Array & vdofs) const; /** S = A_ii - A_ib (A_bb)^{-1} A_bi and y = y_i - A_ib (A_bb)^{-1} y_b. Here (b) and (i) stand for "bubble" and "interface" respectively */ void GetLocalSchurComplement(int el, const Array & tr_idx, const Array & int_idx, const DenseMatrix & elmat, const Vector & elvect, DenseMatrix & rmat, Vector & rvect); void ComputeOffsets(); void BuildProlongation(); #ifdef MFEM_USE_MPI void BuildParallelProlongation(); #endif // ess_tdof list for each space Array *> ess_tdofs; void FillEssTdofLists(const Array & ess_tdof_list); void ConformingAssemble(int skip_zeros); /** Restrict a marker Array on the true FE spaces dofs to a marker Array on the reduced/trace true FE spaces dofs. */ void ConvertMarkerToReducedTrueDofs(Array & tdof_marker, Array & rtdof_marker); void SetSpaces(Array & fes_); void Init(); public: BlockStaticCondensation(Array & fes_); ~BlockStaticCondensation(); /** Assemble the contribution to the Schur complement from the given element matrix 'elmat'; save the other blocks internally: A_bb_inv, A_bi, and A_bi. */ void AssembleReducedSystem(int el, DenseMatrix &elmat, Vector & elvect); /// Finalize the construction of the Schur complement matrix. void Finalize(int skip_zeros = 0); /// Determine and save internally essential reduced true dofs. void SetEssentialTrueDofs(const Array &ess_tdof_list); /// Eliminate the given reduced true dofs from the Schur complement matrix S. void EliminateReducedTrueDofs(const Array &ess_rtdof_list, Matrix::DiagonalPolicy dpolicy); bool HasEliminatedBC() const { #ifndef MFEM_USE_MPI return S_e; #else return S_e || pS_e; #endif } /// Return the serial Schur complement matrix. BlockMatrix &GetSchurMatrix() { return *S; } /// Return the eliminated part of the serial Schur complement matrix. BlockMatrix &GetSchurMatrixElim() { return *S_e; } #ifdef MFEM_USE_MPI /// Return the parallel Schur complement matrix. BlockOperator &GetParallelSchurMatrix() { return *pS; } /// Return the eliminated part of the parallel Schur complement matrix. BlockOperator &GetParallelSchurMatrixElim() { return *pS_e; } void ParallelAssemble(BlockMatrix *m); #endif /** Form the global reduced system matrix using the given @a diag_policy. This method can be called after Assemble() is called. */ void FormSystemMatrix(Operator::DiagonalPolicy diag_policy); /** Restrict a solution vector on the full FE space dofs to a vector on the reduced/trace true FE space dofs. */ void ReduceSolution(const Vector &sol, Vector &sc_sol) const; /** @brief Set the reduced solution `X` and r.h.s `B` vectors from the full linear system solution `x` and r.h.s. `b` vectors. This method should be called after the internal reduced essential dofs have been set using SetEssentialTrueDofs() and both the Schur complement and its eliminated part have been finalized. */ void ReduceSystem(Vector &x, Vector &X, Vector &B, int copy_interior = 0) const; /** Restrict a list of true FE space dofs to a list of reduced/trace true FE space dofs. */ void ConvertListToReducedTrueDofs(const Array &ess_tdof_list, Array &ess_rtdof_list) const; /** Given a solution of the reduced system 'sc_sol' and the RHS 'b' for the full linear system, compute the solution of the full system 'sol'. */ void ComputeSolution(const Vector &sc_sol, Vector &sol) const; }; } #endif