finite_element_vector.hpp

Go to the documentation of this file.

// Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and
// other Serac Project Developers. See the top-level LICENSE file for
// details.
//
// SPDX-License-Identifier: (BSD-3-Clause)
 
#pragma once
 
#include <optional>
 
#include "mfem.hpp"
 
#include "serac/infrastructure/variant.hpp"
#include "serac/numerics/functional/functional.hpp"
 
namespace serac {
 
using GeneralCoefficient = variant<std::shared_ptr<mfem::Coefficient>, std::shared_ptr<mfem::VectorCoefficient>>;
 
enum class ElementType
{
  H1,     
  HCURL,  
  HDIV,   
  L2      
};
 
class FiniteElementVector : public mfem::HypreParVector {
public:
  FiniteElementVector(const mfem::ParFiniteElementSpace& space, const std::string& name = "");
 
  template <typename FunctionSpace>
  FiniteElementVector(mfem::ParMesh& mesh, FunctionSpace, const std::string& name = "") : mesh_(mesh), name_(name)
  {
    const int dim = mesh.Dimension();
 
    const auto ordering = mfem::Ordering::byNODES;
 
    switch (FunctionSpace::family) {
      case Family::H1:
        coll_ = std::make_unique<mfem::H1_FECollection>(FunctionSpace::order, dim);
        break;
      case Family::HCURL:
        coll_ = std::make_unique<mfem::ND_FECollection>(FunctionSpace::order, dim);
        break;
      case Family::HDIV:
        coll_ = std::make_unique<mfem::RT_FECollection>(FunctionSpace::order, dim);
        break;
      case Family::L2:
        // We use GaussLobatto basis functions as this is what is used for the serac::Functional FE kernels
        coll_ = std::make_unique<mfem::L2_FECollection>(FunctionSpace::order, dim, mfem::BasisType::GaussLobatto);
        break;
      default:
        SLIC_ERROR_ROOT("Unknown finite element space requested.");
        break;
    }
 
    space_ = std::make_unique<mfem::ParFiniteElementSpace>(&mesh, coll_.get(), FunctionSpace::components, ordering);
 
    // Construct a hypre par vector based on the new finite element space
    HypreParVector new_vector(space_.get());
 
    // Move the data from this new hypre vector into this object without doubly allocating the data
    auto* parallel_vec = new_vector.StealParVector();
    WrapHypreParVector(parallel_vec);
 
    // Initialize the vector to zero
    HypreParVector::operator=(0.0);
  }
 
  FiniteElementVector(const FiniteElementVector& rhs) : FiniteElementVector(*rhs.space_, rhs.name_)
  {
    HypreParVector::operator=(rhs);
  }
 
  FiniteElementVector(FiniteElementVector&& rhs);
 
  FiniteElementVector& operator=(const FiniteElementVector& rhs);
 
  FiniteElementVector& operator=(FiniteElementVector&& rhs);
 
  FiniteElementVector& operator=(const mfem::HypreParVector& rhs);
 
  FiniteElementVector& operator=(const mfem::Vector& rhs);
 
  MPI_Comm comm() const { return space_->GetComm(); }
 
  mfem::ParMesh& mesh() { return mesh_; }
  const mfem::ParMesh& mesh() const { return mesh_; }
 
  mfem::ParFiniteElementSpace& space() { return *space_; }
  const mfem::ParFiniteElementSpace& space() const { return *space_; }
 
  std::string name() const { return name_; }
 
  FiniteElementVector& operator=(const double value);
 
  virtual ~FiniteElementVector() {}
 
protected:
  std::reference_wrapper<mfem::ParMesh> mesh_;
 
  std::unique_ptr<mfem::FiniteElementCollection> coll_;
 
  std::unique_ptr<mfem::ParFiniteElementSpace> space_;
 
  std::string name_ = "";
};
 
double avg(const FiniteElementVector& fe_vector);
 
double max(const FiniteElementVector& fe_vector);
 
double min(const FiniteElementVector& fe_vector);
 
double innerProduct(const FiniteElementVector& vec1, const FiniteElementVector& vec2);
 
}  // namespace serac