finite_element_state.hpp

Go to the documentation of this file.

// Copyright (c) Lawrence Livermore National Security, LLC and
// other Smith Project Developers. See the top-level LICENSE file for
// details.
//
// SPDX-License-Identifier: (BSD-3-Clause)
 
#pragma once
 
#include <functional>
#include <memory>
#include <optional>
#include <type_traits>
#include <utility>
 
#include "mfem.hpp"
#include "smith/infrastructure/variant.hpp"
#include "smith/physics/state/finite_element_vector.hpp"
#include "smith/numerics/functional/domain.hpp"
#include "smith/numerics/functional/tensor.hpp"
 
namespace smith {
 
namespace detail {
template <int dim>
void setMfemVectorFromTensorOrDouble(mfem::Vector& v_mfem, const tensor<double, dim>& v)
{
  SLIC_ERROR_IF(v_mfem.Size() != dim, "Cannot copy tensor into an MFEM Vector with incompatible size.");
  for (int i = 0; i < dim; i++) v_mfem(i) = v[i];
}
 
inline void setMfemVectorFromTensorOrDouble(mfem::Vector& v_mfem, double v)
{
  SLIC_ERROR_IF(v_mfem.Size() != 1, "Mfem Vector is not a singleton.");
  v_mfem = v;
}
 
template <typename Ret, typename Arg, typename... Rest>
Arg first_argument_helper(Ret (*)(Arg, Rest...));
 
template <typename Ret, typename F, typename Arg, typename... Rest>
Arg first_argument_helper(Ret (F::*)(Arg, Rest...));
 
template <typename Ret, typename F, typename Arg, typename... Rest>
Arg first_argument_helper(Ret (F::*)(Arg, Rest...) const);
 
template <typename F>
decltype(first_argument_helper(&F::operator())) first_argument_helper(F);
 
template <typename T>
using first_argument = std::decay_t<decltype(first_argument_helper(std::declval<T>()))>;
 
template <typename Callable>
auto evaluateTensorFunctionOnMfemVector(const mfem::Vector& X_mfem, Callable&& f)
{
  first_argument<Callable> X;
  SLIC_ERROR_IF(X_mfem.Size() != size(X),
                "Size of tensor in callable does not match spatial dimension of MFEM Vector.");
  for (int i = 0; i < X_mfem.Size(); i++) X[i] = X_mfem[i];
  return f(X);
}
 
}  // namespace detail
 
inline bool is_scalar_valued(const GeneralCoefficient& coef)
{
  return holds_alternative<std::shared_ptr<mfem::Coefficient>>(coef);
}
 
inline bool is_vector_valued(const GeneralCoefficient& coef)
{
  return holds_alternative<std::shared_ptr<mfem::VectorCoefficient>>(coef);
}
 
class FiniteElementState : public FiniteElementVector {
 public:
  using FiniteElementVector::FiniteElementVector;
  using mfem::Vector::Print;
 
  FiniteElementState(const FiniteElementState& rhs) : FiniteElementVector(rhs) {}
 
  FiniteElementState(FiniteElementState&& rhs) : FiniteElementVector(std::move(rhs)) {}
 
  FiniteElementState& operator=(const FiniteElementState& rhs)
  {
    FiniteElementVector::operator=(rhs);
    return *this;
  }
 
  FiniteElementState& operator=(FiniteElementState&& rhs)
  {
    FiniteElementVector::operator=(rhs);
    return *this;
  }
 
  FiniteElementState& operator=(const mfem::HypreParVector& rhs)
  {
    FiniteElementVector::operator=(rhs);
    return *this;
  }
 
  FiniteElementState& operator=(const mfem::Vector& rhs)
  {
    FiniteElementVector::operator=(rhs);
    return *this;
  }
 
  FiniteElementState& operator=(double rhs)
  {
    FiniteElementVector::operator=(rhs);
    return *this;
  }
 
  void fillGridFunction(mfem::ParGridFunction& grid_function) const { grid_function.SetFromTrueDofs(*this); }
 
  void setFromGridFunction(const mfem::ParGridFunction& grid_function) { grid_function.GetTrueDofs(*this); }
 
  void project(mfem::VectorCoefficient& coef, mfem::Array<int>& dof_list);
 
  void project(mfem::Coefficient& coef, mfem::Array<int>& dof_list, std::optional<int> component = {});
 
  void project(const GeneralCoefficient& coef);
 
  void project(mfem::Coefficient& coef);
 
  void project(mfem::VectorCoefficient& coef);
 
  void projectOnBoundary(mfem::Coefficient& coef, const mfem::Array<int>& markers);
 
  void projectOnBoundary(mfem::VectorCoefficient& coef, const mfem::Array<int>& markers);
 
  void project(mfem::Coefficient& coef, const Domain& d);
 
  void project(mfem::VectorCoefficient& coef, const Domain& d);
 
  template <typename FieldFunction>
  void setFromFieldFunction(FieldFunction&& field_function)
  {
    auto evaluate_mfem = [&field_function](const mfem::Vector& X_mfem, mfem::Vector& u_mfem) {
      auto u = detail::evaluateTensorFunctionOnMfemVector(X_mfem, field_function);
      detail::setMfemVectorFromTensorOrDouble(u_mfem, u);
    };
 
    mfem::VectorFunctionCoefficient coef(space_->GetVDim(), evaluate_mfem);
    project(coef);
  }
 
  mfem::ParGridFunction& gridFunction() const;
 
 protected:
  mutable std::unique_ptr<mfem::ParGridFunction> grid_func_;
};
 
double norm(const FiniteElementState& state, const double p = 2);
 
double computeL2Error(const FiniteElementState& state, mfem::VectorCoefficient& exact_solution);
 
double computeL2Error(const FiniteElementState& state, mfem::Coefficient& exact_solution);
 
}  // namespace smith