/**
 * \file
 * \copyright
 * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
 *            Distributed under a Modified BSD License.
 *              See accompanying file LICENSE.txt or
 *              http://www.opengeosys.org/project/license
 *
 */

#pragma once

#include <memory>
#include <vector>
#include "MaterialLib/Fluid/FluidPropertyHeaders.h"
#include "MaterialLib/PhysicalConstant.h"
#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
#include "MaterialLib/PorousMedium/PorousPropertyHeaders.h"
#include "MaterialLib/PorousMedium/Storage/Storage.h"
#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h"
#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.h"
#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.h"
#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/RelativePermeability.h"
#include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"

namespace MeshLib
{
template <typename PROP_VAL_TYPE>
class PropertyVector;
}

namespace ProcessLib::TwoPhaseFlowWithPrho
{
class TwoPhaseFlowWithPrhoMaterialProperties
{
public:
    using ArrayType = MaterialLib::Fluid::FluidProperty::ArrayType;

    TwoPhaseFlowWithPrhoMaterialProperties(
        boost::optional<MeshLib::PropertyVector<int> const&> const material_ids,
        std::unique_ptr<MaterialLib::Fluid::FluidProperty>
            liquid_density,
        std::unique_ptr<MaterialLib::Fluid::FluidProperty>
            viscosity,
        std::unique_ptr<MaterialLib::Fluid::FluidProperty>
            gas_density,
        std::unique_ptr<MaterialLib::Fluid::FluidProperty>
            gas_viscosity,
        std::vector<std::unique_ptr<MaterialLib::PorousMedium::Permeability>>&&
            intrinsic_permeability_models,
        std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>>&&
            porosity_models,
        std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>>&&
            storage_models,
        std::vector<std::unique_ptr<
            MaterialLib::PorousMedium::CapillaryPressureSaturation>>&&
            capillary_pressure_models,
        std::vector<
            std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>>&&
            relative_permeability_models);

    int getMaterialID(const std::size_t element_id);

    Eigen::MatrixXd getPermeability(const int material_id,
                                    const double t,
                                    const ParameterLib::SpatialPosition& pos,
                                    const int dim) const;

    double getPorosity(const int material_id, const double t,
                       const ParameterLib::SpatialPosition& pos, const double p,
                       const double T, const double porosity_variable) const;

    double getNonwetRelativePermeability(
        const double t, const ParameterLib::SpatialPosition& pos,
        const double p, const double T, const double saturation) const;
    double getWetRelativePermeability(const double t,
                                      const ParameterLib::SpatialPosition& pos,
                                      const double p, const double T,
                                      const double saturation) const;
    double getCapillaryPressure(const int material_id, const double t,
                                const ParameterLib::SpatialPosition& pos,
                                const double p, const double T,
                                const double saturation) const;
    double getCapillaryPressureDerivative(
        const int material_id, const double t,
        const ParameterLib::SpatialPosition& pos, const double p,
        const double T, const double saturation) const;
    double getLiquidDensity(const double p, const double T) const;
    double getGasDensity(const double p, const double T) const;
    double getGasViscosity(const double p, const double T) const;
    double getLiquidViscosity(const double p, const double T) const;
    bool computeConstitutiveRelation(
        double const t,
        ParameterLib::SpatialPosition const& x_position,
        const int material_id,
        double const pg,
        double const X,
        double const T,
        double& Sw,
        double& X_m,
        double& dsw_dpg,
        double& dsw_dX,
        double& dxm_dpg,
        double& dxm_dX);

protected:
    std::unique_ptr<MaterialLib::Fluid::FluidProperty> _liquid_density;
    std::unique_ptr<MaterialLib::Fluid::FluidProperty> _viscosity;
    std::unique_ptr<MaterialLib::Fluid::FluidProperty> _gas_density;
    std::unique_ptr<MaterialLib::Fluid::FluidProperty> _gas_viscosity;

    /** Use two phase models for different material zones.
    *  Material IDs must be given as mesh element properties.
    */
    boost::optional<MeshLib::PropertyVector<int> const&> const _material_ids;

    std::vector<std::unique_ptr<MaterialLib::PorousMedium::Permeability>>
        _intrinsic_permeability_models;
    std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>>
        _porosity_models;
    std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>>
        _storage_models;
    std::vector<
        std::unique_ptr<MaterialLib::PorousMedium::CapillaryPressureSaturation>>
        _capillary_pressure_models;
    std::vector<
        std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>>
        _relative_permeability_models;

private:
    static int const jacobian_residual_size = 2;
    using ResidualVector = Eigen::Matrix<double, jacobian_residual_size, 1>;
    using JacobianMatrix =
        Eigen::Matrix<double, jacobian_residual_size, jacobian_residual_size,
                      Eigen::RowMajor>;
    using UnknownVector = Eigen::Matrix<double, jacobian_residual_size, 1>;

private:
    /**
    * Calculates the residual vector.
    */
    void calculateResidual(const int material_id, double const pl,
                           double const X, double const T, double Sw,
                           double rho_h2_wet, ResidualVector& res);
    /**
    * Calculates the Jacobian.
    */
    void calculateJacobian(const int material_id, double const t,
                           ParameterLib::SpatialPosition const& x,
                           double const pl, double const X, double const T,
                           JacobianMatrix& Jac, double Sw, double rho_h2_wet);
    /** Complementary condition 1
    * for calculating molar fraction of light component in the liquid phase
    */
    static double calculateEquilibiumRhoWetLight(double const pg,
                                                 double const Sw,
                                                 double const rho_wet_h2);
    /** Complementary condition 2
    * for calculating the saturation
    */
    static double calculateSaturation(double /*PL*/, double X, double Sw,
                                      double rho_wet_h2, double rho_nonwet_h2,
                                      double /*T*/);
    /**
    * Calculate the derivatives using the analytical way
    */
    double calculatedSwdP(double pl, double S, double rho_wet_h2,
                          double const T, int current_material_id) const;
    /**
    * Calculate the derivatives using the analytical way
    */
    double calculatedSwdX(double const pl, const double /*X*/, const double S,
                          const double rho_wet_h2, double const T,
                          int current_material_id) const;
    /**
    * Calculate the derivatives using the analytical way
    */
    double calculatedXmdX(double pl, double Sw, double rho_wet_h2, double dSwdX,
                          int current_material_id) const;
    /**
    * Calculate the derivatives using the analytical way
    */
    double calculatedXmdP(double pl, double Sw, double rho_wet_h2, double dSwdP,
                          int current_material_id) const;
};

}  // namespace ProcessLib::TwoPhaseFlowWithPrho