// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program 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 General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
/*!
 * \file
 *
 * \brief The spatial parameters for the exercise-fluidsystem problem
 *        which uses the two-phase and two-phase two-component box model.
 */
#ifndef DUMUX_EXERCISE_FLUIDSYSTEM_SPATIAL_PARAMS_HH
#define DUMUX_EXERCISE_FLUIDSYSTEM_SPATIAL_PARAMS_HH

// include parent spatialparameters
#include <dumux/material/spatialparams/fv.hh>

//inlclude fluid matrix interaction relationship
#include <dumux/material/fluidmatrixinteractions/2p/brookscorey.hh>
#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>

namespace Dumux {

/*!
 * \ingroup TwoPBoxModel
 *
 * \brief The spatial parameters for the exercise-fluidsystem problem
 *        which uses the two-phase and two-phase two-component box model.
 */
template<class FVGridGeometry, class Scalar>
class ExerciseFluidsystemSpatialParams
: public FVSpatialParams<FVGridGeometry, Scalar, ExerciseFluidsystemSpatialParams<FVGridGeometry, Scalar>>
{
    using ThisType = ExerciseFluidsystemSpatialParams<FVGridGeometry, Scalar>;
    using ParentType = FVSpatialParams<FVGridGeometry, Scalar, ThisType>;
    using GridView = typename FVGridGeometry::GridView;

    static constexpr int dim = GridView::dimension;
    static constexpr int dimWorld = GridView::dimensionworld;
    using Element = typename GridView::template Codim<0>::Entity;
    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;

    using PcKrSwCurve = FluidMatrix::BrooksCoreyDefault<Scalar>;

public:
    // export permeability type
    using PermeabilityType = Dune::FieldMatrix<Scalar, dim, dim>;

    /*!
     * \brief The constructor
     *
     * \param fvGridGeometry The finite volume grid geometry
     */
    ExerciseFluidsystemSpatialParams(std::shared_ptr<const FVGridGeometry>& fvGridGeometry)
    : ParentType(fvGridGeometry)
    , K_(0)
    , KLens_(0)
    , pcKrSwCurve_("SpatialParams")
    , lensPcKrSwCurve_("SpatialParams.Lens")
    {
        //set main diagonal entries of the permeability tensor to a value
        //setting to one value means: isotropic, homogeneous
        for (int i = 0; i < dim; i++)
        {
            K_[i][i] = 1e-7;
            KLens_[i][i] = 1e-10;
        }
    }


    /*!
     * \brief Define the intrinsic permeability \f$\mathrm{[m^2]}\f$.
     *
     * \param globalPos The global position
     */
    PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const

    {
        if (isInLens(globalPos))
            return KLens_;
        return K_;
    }

    /*!
     * \brief Define the porosity \f$\mathrm{[-]}\f$.
     *
     * \param globalPos The global position
     */
    Scalar porosityAtPos(const GlobalPosition& globalPos) const
    {
        if (isInLens(globalPos))
            return 0.1;
        return 0.2;
    }

    /*!
     * \brief Returns the fluid-matrix interaction law at a given location
     * \param globalPos The global coordinates for the given location
     */
    auto fluidMatrixInteractionAtPos(const GlobalPosition& globalPos) const
    {
        if (isInLens(globalPos))
            return makeFluidMatrixInteraction(pcKrSwCurve_);
        return makeFluidMatrixInteraction(lensPcKrSwCurve_);
    }

    /*!
     * \brief Function for defining which phase is to be considered as the wetting phase.
     *
     * \return the wetting phase index
     * \param globalPos The position of the center of the element
     */
    template<class FluidSystem>
    int wettingPhaseAtPos(const GlobalPosition& globalPos) const
    {
        // Our fluid system is H2OMyCompressibleComponent
        // We want to define water as the wetting phase in
        // the entire domain (see fluid system for the phase indices)

        // TODO: dumux-course-task 4
        // Adapt the following line so that the phase of our new component is
        // the wetting phase, only within the lenses.
        return FluidSystem::phase0Idx;
    }

    //! if we are in the lens
    bool isInLens(const GlobalPosition& globalPos) const
    {
        const auto x = globalPos[0];
        const auto y = globalPos[1];
        return (x < 40 + eps_ && x > 20 - eps_ && y > 35 - eps_ && y < 45 + eps_) ||
               (x < 50 + eps_ && x > 30 - eps_ && y < 30 + eps_ && y > 15 - eps_);
    }

private:
    static constexpr Scalar eps_ = 1e-6;

    Dune::FieldMatrix<Scalar, dim, dim> K_;
    Dune::FieldMatrix<Scalar, dim, dim> KLens_;

    const PcKrSwCurve pcKrSwCurve_;
    const PcKrSwCurve lensPcKrSwCurve_;
};
} // end namespace Dumux
#endif