constraint.hpp

Go to the documentation of this file.

// Copyright 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 <vector>
#include "smith/physics/common.hpp"
#include "smith/physics/field_types.hpp"
 
namespace mfem {
class Vector;
class HypreParMatrix;
}  // namespace mfem
 
namespace smith {
 
class FiniteElementState;
 
class Constraint {
 public:
  Constraint(const std::string& name = "constraint") : name_(name) {}
 
  virtual ~Constraint() {}
 
  virtual mfem::Vector evaluate(double time, double dt, const std::vector<ConstFieldPtr>& fields) const = 0;
 
  virtual std::unique_ptr<mfem::HypreParMatrix> jacobian(double time, double dt,
                                                         const std::vector<ConstFieldPtr>& fields,
                                                         int direction) const = 0;
 
  virtual std::unique_ptr<mfem::HypreParMatrix> jacobian_tilde(double time, double dt,
                                                               const std::vector<ConstFieldPtr>& fields,
                                                               int direction) const
  {
    return jacobian(time, dt, fields, direction);
  };
 
  virtual mfem::Vector residual_contribution(double time, double dt, const std::vector<ConstFieldPtr>& fields,
                                             const mfem::Vector& multipliers, int direction) const
  {
    std::unique_ptr<mfem::HypreParMatrix> jac = jacobian_tilde(time, dt, fields, direction);
    mfem::Vector y(jac->Width());
    y = 0.0;
    SLIC_ERROR_ROOT_IF(jac->Height() != multipliers.Size(), "Incompatible matrix and vector sizes.");
    jac->MultTranspose(multipliers, y);
    return y;
  };
 
  virtual std::unique_ptr<mfem::HypreParMatrix> residual_contribution_jacobian(
      [[maybe_unused]] double time, [[maybe_unused]] double dt,
      [[maybe_unused]] const std::vector<ConstFieldPtr>& fields, [[maybe_unused]] const mfem::Vector& multipliers,
      [[maybe_unused]] int direction) const
  {
    SLIC_ERROR_ROOT(axom::fmt::format("Base class must override residual_contribution_jacobian before usage"));
    std::unique_ptr<mfem::HypreParMatrix> res_contr_jacobian = nullptr;
    return res_contr_jacobian;
  };
 
  std::string name() const { return name_; }
 
 private:
  std::string name_;
};
 
}  // namespace smith