Line data    Source code

       1              : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
       2              : // SPDX-License-Identifier: Apache-2.0
       3              : 
       4              : #ifndef PALACE_LINALG_DIV_FREE_HPP
       5              : #define PALACE_LINALG_DIV_FREE_HPP
       6              : 
       7              : #include <memory>
       8              : #include <vector>
       9              : #include "linalg/ksp.hpp"
      10              : #include "linalg/operator.hpp"
      11              : #include "linalg/vector.hpp"
      12              : 
      13              : namespace mfem
      14              : {
      15              : 
      16              : template <typename T>
      17              : class Array;
      18              : 
      19              : }  // namespace mfem
      20              : 
      21              : namespace palace
      22              : {
      23              : 
      24              : class FiniteElementSpaceHierarchy;
      25              : class FiniteElementSpace;
      26              : class MaterialOperator;
      27              : 
      28              : //
      29              : // This solver implements a projection onto a divergence-free space satisfying Gᵀ M x = 0,
      30              : // where G represents the discrete gradient matrix with columns spanning the nullspace of
      31              : // the curl-curl operator.
      32              : //
      33              : template <typename VecType>
      34              : class DivFreeSolver
      35              : {
      36              :   using OperType = typename std::conditional<std::is_same<VecType, ComplexVector>::value,
      37              :                                              ComplexOperator, Operator>::type;
      38              : 
      39              : private:
      40              :   // Operators for the divergence-free projection.
      41              :   std::unique_ptr<OperType> M;
      42              :   std::unique_ptr<Operator> WeakDiv;
      43              :   const Operator *Grad;
      44              :   const mfem::Array<int> *bdr_tdof_list_M;
      45              : 
      46              :   // Optional storage for homogeneous Dirichlet boundary condition on a single true dof,
      47              :   // used when the input array of H1 boundary dofs is empty to prevent the Poisson operator
      48              :   // from being singular.
      49              :   std::vector<mfem::Array<int>> aux_tdof_lists;
      50              : 
      51              :   // Linear solver for the projected linear system (Gᵀ M G) y = x.
      52              :   std::unique_ptr<BaseKspSolver<OperType>> ksp;
      53              : 
      54              :   // Workspace objects for solver application.
      55              :   mutable VecType psi, rhs;
      56              : 
      57              : public:
      58              :   DivFreeSolver(const MaterialOperator &mat_op, FiniteElementSpace &nd_fespace,
      59              :                 FiniteElementSpaceHierarchy &h1_fespaces,
      60              :                 const std::vector<mfem::Array<int>> &h1_bdr_tdof_lists, double tol,
      61              :                 int max_it, int print);
      62              : 
      63              :   // Given a vector of Nedelec dofs for an arbitrary vector field, compute the Nedelec dofs
      64              :   // of the irrotational portion of this vector field. The resulting vector will satisfy
      65              :   // ∇ x y = 0.
      66              :   void Mult(VecType &y) const;
      67              : 
      68            0 :   void Mult(const VecType &x, VecType &y) const
      69              :   {
      70            0 :     y = x;
      71            0 :     Mult(y);
      72            0 :   }
      73              : };
      74              : 
      75              : }  // namespace palace
      76              : 
      77              : #endif  // PALACE_LINALG_DIV_FREE_HPP