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_MULTIGRID_HPP
       5              : #define PALACE_FEM_MULTIGRID_HPP
       6              : 
       7              : #include <memory>
       8              : #include <vector>
       9              : #include <mfem.hpp>
      10              : #include "fem/fespace.hpp"
      11              : #include "fem/mesh.hpp"
      12              : #include "utils/geodata.hpp"
      13              : #include "utils/iodata.hpp"
      14              : 
      15              : namespace palace::fem
      16              : {
      17              : 
      18              : //
      19              : // Methods for constructing hierarchies of finite element spaces for geometric multigrid.
      20              : //
      21              : 
      22              : // Construct sequence of FECollection objects.
      23              : template <typename FECollection>
      24              : inline std::vector<std::unique_ptr<FECollection>>
      25           63 : ConstructFECollections(int p, int dim, int mg_max_levels, MultigridCoarsening mg_coarsening,
      26              :                        bool mat_lor)
      27              : {
      28              :   // If the solver will use a LOR preconditioner, we need to construct with a specific basis
      29              :   // type.
      30              :   constexpr int pmin = (std::is_base_of<mfem::H1_FECollection, FECollection>::value ||
      31              :                         std::is_base_of<mfem::ND_FECollection, FECollection>::value)
      32              :                            ? 1
      33              :                            : 0;
      34           63 :   MFEM_VERIFY(p >= pmin, "FE space order must not be less than " << pmin << "!");
      35           63 :   int b1 = mfem::BasisType::GaussLobatto, b2 = mfem::BasisType::GaussLegendre;
      36           40 :   if (mat_lor)
      37              :   {
      38            0 :     b2 = mfem::BasisType::IntegratedGLL;
      39              :   }
      40              : 
      41              :   // Construct the p-multigrid hierarchy, first finest to coarsest and then reverse the
      42              :   // order.
      43              :   std::vector<std::unique_ptr<FECollection>> fecs;
      44           63 :   for (int l = 0; l < std::max(1, mg_max_levels); l++)
      45              :   {
      46              :     if constexpr (std::is_base_of<mfem::ND_FECollection, FECollection>::value ||
      47              :                   std::is_base_of<mfem::RT_FECollection, FECollection>::value)
      48              :     {
      49           40 :       fecs.push_back(std::make_unique<FECollection>(p, dim, b1, b2));
      50              :     }
      51              :     else
      52              :     {
      53           23 :       fecs.push_back(std::make_unique<FECollection>(p, dim, b1));
      54              :       MFEM_CONTRACT_VAR(b2);
      55              :     }
      56           63 :     if (p == pmin)
      57              :     {
      58              :       break;
      59              :     }
      60            0 :     switch (mg_coarsening)
      61              :     {
      62            0 :       case MultigridCoarsening::LINEAR:
      63            0 :         p--;
      64            0 :         break;
      65            0 :       case MultigridCoarsening::LOGARITHMIC:
      66            0 :         p = (p + pmin) / 2;
      67            0 :         break;
      68              :     }
      69              :   }
      70              :   std::reverse(fecs.begin(), fecs.end());
      71              : 
      72           63 :   return fecs;
      73            0 : }
      74              : 
      75              : // Construct a hierarchy of finite element spaces given a sequence of meshes and finite
      76              : // element collections. Additionally, Dirichlet boundary conditions are marked.
      77              : template <typename FECollection>
      78           63 : inline FiniteElementSpaceHierarchy ConstructFiniteElementSpaceHierarchy(
      79              :     int mg_max_levels, const std::vector<std::unique_ptr<Mesh>> &mesh,
      80              :     const std::vector<std::unique_ptr<FECollection>> &fecs,
      81              :     const mfem::Array<int> *dbc_attr = nullptr,
      82              :     std::vector<mfem::Array<int>> *dbc_tdof_lists = nullptr)
      83              : {
      84           63 :   MFEM_VERIFY(!mesh.empty() && !fecs.empty() &&
      85              :                   (!dbc_tdof_lists || dbc_tdof_lists->empty()),
      86              :               "Empty mesh or FE collection for FE space construction!");
      87           63 :   int coarse_mesh_l = std::max(0, static_cast<int>(mesh.size() + fecs.size()) - 1 -
      88              :                                       std::max(1, mg_max_levels));
      89           63 :   FiniteElementSpaceHierarchy fespaces(
      90          126 :       std::make_unique<FiniteElementSpace>(*mesh[coarse_mesh_l], fecs[0].get()));
      91              : 
      92              :   mfem::Array<int> dbc_marker;
      93           63 :   if (dbc_attr && dbc_tdof_lists)
      94              :   {
      95              :     int bdr_attr_max = mesh[coarse_mesh_l]->Get().bdr_attributes.Size()
      96           40 :                            ? mesh[coarse_mesh_l]->Get().bdr_attributes.Max()
      97              :                            : 0;
      98           80 :     dbc_marker = mesh::AttrToMarker(bdr_attr_max, *dbc_attr);
      99           40 :     fespaces.GetFinestFESpace().Get().GetEssentialTrueDofs(dbc_marker,
     100           40 :                                                            dbc_tdof_lists->emplace_back());
     101              :   }
     102              : 
     103              :   // h-refinement.
     104           63 :   for (std::size_t l = coarse_mesh_l + 1; l < mesh.size(); l++)
     105              :   {
     106            0 :     fespaces.AddLevel(std::make_unique<FiniteElementSpace>(*mesh[l], fecs[0].get()));
     107            0 :     if (dbc_attr && dbc_tdof_lists)
     108              :     {
     109            0 :       fespaces.GetFinestFESpace().Get().GetEssentialTrueDofs(
     110            0 :           dbc_marker, dbc_tdof_lists->emplace_back());
     111              :     }
     112              :   }
     113              : 
     114              :   // p-refinement.
     115           63 :   for (std::size_t l = 1; l < fecs.size(); l++)
     116              :   {
     117            0 :     fespaces.AddLevel(std::make_unique<FiniteElementSpace>(*mesh.back(), fecs[l].get()));
     118            0 :     if (dbc_attr && dbc_tdof_lists)
     119              :     {
     120            0 :       fespaces.GetFinestFESpace().Get().GetEssentialTrueDofs(
     121            0 :           dbc_marker, dbc_tdof_lists->emplace_back());
     122              :     }
     123              :   }
     124              : 
     125           63 :   return fespaces;
     126              : }
     127              : 
     128              : }  // namespace palace::fem
     129              : 
     130              : #endif  // PALACE_FEM_MULTIGRID_HPP