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_FEM_BILINEARFORM_HPP
       5              : #define PALACE_FEM_BILINEARFORM_HPP
       6              : 
       7              : #include <memory>
       8              : #include <vector>
       9              : #include <mfem.hpp>
      10              : #include "fem/integrator.hpp"
      11              : #include "fem/libceed/operator.hpp"
      12              : #include "linalg/hypre.hpp"
      13              : 
      14              : namespace palace
      15              : {
      16              : 
      17              : class FiniteElementSpace;
      18              : class FiniteElementSpaceHierarchy;
      19              : 
      20              : //
      21              : // This class implements bilinear and mixed bilinear forms based on integrators assembled
      22              : // using the libCEED library. Assembly in the form of a partially assembled operator or
      23              : // fully assembled sparse matrix is available.
      24              : //
      25        11100 : class BilinearForm
      26              : {
      27              : protected:
      28              :   // Domain and range finite element spaces.
      29              :   const FiniteElementSpace &trial_fespace, &test_fespace;
      30              : 
      31              :   // List of domain and boundary integrators making up the bilinear form.
      32              :   std::vector<std::unique_ptr<BilinearFormIntegrator>> domain_integs, boundary_integs;
      33              : 
      34              :   std::unique_ptr<ceed::Operator>
      35              :   PartialAssemble(const FiniteElementSpace &trial_fespace,
      36              :                   const FiniteElementSpace &test_fespace) const;
      37              : 
      38              : public:
      39              :   // Order above which to use partial assembly vs. full.
      40              :   inline static int pa_order_threshold = 1;
      41              : 
      42              : public:
      43              :   BilinearForm(const FiniteElementSpace &trial_fespace,
      44              :                const FiniteElementSpace &test_fespace)
      45        11100 :     : trial_fespace(trial_fespace), test_fespace(test_fespace)
      46              :   {
      47              :   }
      48              :   BilinearForm(const FiniteElementSpace &fespace) : BilinearForm(fespace, fespace) {}
      49              : 
      50        27294 :   const auto &GetTrialSpace() const { return trial_fespace; }
      51        16548 :   const auto &GetTestSpace() const { return test_fespace; }
      52              : 
      53              :   template <typename T, typename... U>
      54         4596 :   void AddDomainIntegrator(U &&...args)
      55              :   {
      56         4596 :     domain_integs.push_back(std::make_unique<T>(std::forward<U>(args)...));
      57         4596 :   }
      58              : 
      59              :   template <typename T, typename... U>
      60         1644 :   void AddBoundaryIntegrator(U &&...args)
      61              :   {
      62         1644 :     boundary_integs.push_back(std::make_unique<T>(std::forward<U>(args)...));
      63         1644 :   }
      64              : 
      65              :   void AssembleQuadratureData();
      66              : 
      67              :   std::unique_ptr<ceed::Operator> PartialAssemble() const
      68              :   {
      69        10788 :     return PartialAssemble(GetTrialSpace(), GetTestSpace());
      70              :   }
      71              : 
      72            0 :   std::unique_ptr<hypre::HypreCSRMatrix> FullAssemble(bool skip_zeros) const
      73              :   {
      74            0 :     return FullAssemble(*PartialAssemble(), skip_zeros, false);
      75              :   }
      76              : 
      77              :   static std::unique_ptr<hypre::HypreCSRMatrix> FullAssemble(const ceed::Operator &op,
      78              :                                                              bool skip_zeros)
      79              :   {
      80        10710 :     return FullAssemble(op, skip_zeros, false);
      81              :   }
      82              : 
      83              :   static std::unique_ptr<hypre::HypreCSRMatrix> FullAssemble(const ceed::Operator &op,
      84              :                                                              bool skip_zeros, bool set);
      85              : 
      86              :   std::unique_ptr<Operator> Assemble(bool skip_zeros) const;
      87              : 
      88              :   std::vector<std::unique_ptr<Operator>>
      89              :   Assemble(const FiniteElementSpaceHierarchy &fespaces, bool skip_zeros,
      90              :            std::size_t l0 = 0) const;
      91              : };
      92              : 
      93              : // Discrete linear operators map primal vectors to primal vectors for interpolation between
      94              : // spaces.
      95          558 : class DiscreteLinearOperator
      96              : {
      97              : private:
      98              :   // Domain and range finite element spaces.
      99              :   const FiniteElementSpace &trial_fespace, &test_fespace;
     100              : 
     101              :   // List of domain interpolators making up the discrete linear operator.
     102              :   std::vector<std::unique_ptr<DiscreteInterpolator>> domain_interps;
     103              : 
     104              : public:
     105              :   DiscreteLinearOperator(const FiniteElementSpace &trial_fespace,
     106              :                          const FiniteElementSpace &test_fespace)
     107          558 :     : trial_fespace(trial_fespace), test_fespace(test_fespace)
     108              :   {
     109              :   }
     110              : 
     111          216 :   const auto &GetTrialSpace() const { return trial_fespace; }
     112          216 :   const auto &GetTestSpace() const { return test_fespace; }
     113              : 
     114              :   template <typename T, typename... U>
     115          504 :   void AddDomainInterpolator(U &&...args)
     116              :   {
     117          504 :     domain_interps.push_back(std::make_unique<T>(std::forward<U>(args)...));
     118          504 :   }
     119              : 
     120              :   std::unique_ptr<ceed::Operator> PartialAssemble() const;
     121              : 
     122              :   std::unique_ptr<hypre::HypreCSRMatrix> FullAssemble(bool skip_zeros) const
     123              :   {
     124              :     return BilinearForm::FullAssemble(*PartialAssemble(), skip_zeros, true);
     125              :   }
     126              : 
     127              :   static std::unique_ptr<hypre::HypreCSRMatrix> FullAssemble(const ceed::Operator &op,
     128              :                                                              bool skip_zeros)
     129              :   {
     130          222 :     return BilinearForm::FullAssemble(op, skip_zeros, true);
     131              :   }
     132              : };
     133              : 
     134              : }  // namespace palace
     135              : 
     136              : #endif  // PALACE_FEM_BILINEARFORM_HPP