integrators.h

Go to the documentation of this file.

/* ***********************************************************************
//
//   Copyright (C) 2025 -- The 4D-STAR Collaboration
//   File Author: Emily Boudreaux
//   Last Modified: April 21, 2025
//
//   4DSSE is free software; you can use it and/or modify
//   it under the terms and restrictions the GNU General Library Public
//   License version 3 (GPLv3) as published by the Free Software Foundation.
//
//   4DSSE is distributed in the hope that it will be useful,
//   but WITHOUT ANY WARRANTY; without even the implied warranty of
//   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
//   See the GNU Library General Public License for more details.
//
//   You should have received a copy of the GNU Library General Public License
//   along with this software; if not, write to the Free Software
//   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#pragma once
 
#include "mfem.hpp"
#include <string>
#include "config.h"
#include "probe.h"
 

 
namespace serif::polytrope {
  namespace polyMFEMUtils {
    class NonlinearPowerIntegrator: public mfem::NonlinearFormIntegrator {
    public:
      NonlinearPowerIntegrator(double n);

      virtual void AssembleElementVector(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::Vector &elvect) override;
      virtual void AssembleElementGrad (const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::DenseMatrix &elmat) override;
    private:
      serif::config::Config& m_config = serif::config::Config::getInstance();
      serif::probe::LogManager& m_logManager = serif::probe::LogManager::getInstance();
      quill::Logger* m_logger = m_logManager.getLogger("log");
      double m_polytropicIndex;
      double m_epsilon;
      static constexpr double m_regularizationRadius = 0.15; 
      static constexpr double m_regularizationCoeff = 1.0/6.0; 
    };
 
    inline double dfmod(const double epsilon, const double n) {
      if (n == 0.0) {
        return 0.0;
      }
      if (n == 1.0) {
        return 1.0;
      }
      return n * std::pow(epsilon, n - 1.0);
    }
 
    inline double fmod(const double theta, const double n, const double epsilon) {
      if (n == 0.0) {
        return 1.0;
      }
      // For n != 0
      const double y_prime_at_epsilon = dfmod(epsilon, n); // Uses the robust dfmod
      const double y_at_epsilon = std::pow(epsilon, n);     // epsilon^n
 
      // f_mod(theta) = y_at_epsilon + y_prime_at_epsilon * (theta - epsilon)
      return y_at_epsilon + y_prime_at_epsilon * (theta - epsilon);
    }
 
 
  } // namespace polyMFEMUtils
}