#include "MeshParameters.h"
#include <Bnd_Box.hxx>
#include <BRepBndLib.hxx>
#include <TopExp_Explorer.hxx>
#include <GProp_GProps.hxx>
#include <BRepGProp.hxx>
#include <algorithm>

MeshParameters::MeshParameters(ShapeType st)
{
    density_factor = 1.0;
    maxh = 1000;
    minh = 0;
    grading = 0.4;
    param_type = st;
}

MeshParameters::~MeshParameters()
{
    //delete meshing_parameters;
}

MeshParameters::MeshParameters(TopoDS_Shape partition)
{
    //calculate minimal and maximal distances
    Bnd_Box box;
    box.SetGap(0.0);
    BRepBndLib::Add(partition, box);

    Standard_Real xmin, ymin, zmin, xmax, ymax, zmax;
    box.Get(xmin, ymin, zmin, xmax, ymax, zmax);

    double maxdist = std::sqrt(pow((xmax - xmin), 2) + pow((ymax - ymin), 2) + pow((zmax - zmin), 2));

    double mindist = maxdist;

    for (TopExp_Explorer edgeExplorer(partition, TopAbs_ShapeEnum::TopAbs_EDGE); edgeExplorer.More(); edgeExplorer.Next())
    {
        TopoDS_Shape thisEdge = edgeExplorer.Current();
        GProp_GProps massProps;
        BRepGProp::LinearProperties(thisEdge, massProps);
        if (massProps.Mass() < mindist)
            mindist = massProps.Mass();
    }

    density_factor = 1.0;
    minh = std::max(mindist, maxdist / 500);
    maxh = maxdist / 20;
    grading = 0.4;
}

MeshParameters::MeshParameters(const MeshParameters& that)
{
    density_factor = that.density_factor;
    minh = that.minh;
    maxh = that.maxh;
    grading = that.grading;
    param_type = that.param_type;
}

void MeshParameters::operator = (const MeshParameters& that)
{
    density_factor = that.density_factor;
    minh = that.minh;
    maxh = that.maxh;
    grading = that.grading;
    param_type = that.param_type;
}

void MeshParameters::setMaxH(double s)
{
    maxh = s;
}

void MeshParameters::setMinH(double s)
{
    minh = s;
}

void MeshParameters::setGrading(double s)
{
    grading = s;
}

void MeshParameters::setDensityFactor(double s)
{
    density_factor = s;
}

double MeshParameters::getMaxH()
{
    return maxh;
}
double MeshParameters::getGrading()
{
    return grading;
}

double MeshParameters::getMinH()
{
    return minh;
}

double MeshParameters::getDensityFactor()
{
    return density_factor;
}

void MeshParameters::calculateLayers(double skin_min, double skin_max, double geom_limit)
{
    double skin_factor = 1.55; // multiply max skin depth by this number (empirical result)
    double layer_factor = 1.4; // layer growth factor

    // check if max skin depth is larger than geometrical limit
    if ( (geom_limit != 0) && (geom_limit < skin_max * skin_factor) )
    {
        //reset skin_max to geometrical limit
        //this is done to restrict possibility that skin size is larger than thickness of object
        skin_max = geom_limit;
        skin_factor = 1.0;
    }

    //start with 7 layers
    int nr_layers = 7;
    double coeff_sum = 0;
    for (int i = 0; i < nr_layers; i++)
        coeff_sum = coeff_sum + std::pow(layer_factor, i);

    double first_layer = skin_max * skin_factor / coeff_sum;

    // coefficient for 4 layers
    double coeff_4 = 0;
    for (int i = 0; i < 4; i++)
        coeff_4 = coeff_4 + std::pow(layer_factor, i);

    // increase number of layers while minimal thickness has less than 4 layers
    while (skin_min < first_layer * coeff_4)
    {
        nr_layers++;

        //recalculate first layer thickness
        coeff_sum = 0;
        for (int i = 0; i < nr_layers; i++)
            coeff_sum = coeff_sum + std::pow(layer_factor, i);

        first_layer = skin_max * skin_factor / coeff_sum;

        // Should this limitation be here?
        if (nr_layers >= 20)
            break;
    }

    thicknesses = {};;
    for (int i = 0; i < nr_layers; i++)
        thicknesses.push_back(first_layer * std::pow(layer_factor, i));
}

SalomeVLParams MeshParameters::getSalomeVL(double skin_min, double skin_max)
{
    double skin_factor = 1.55; // multiply max skin depth by this number (empirical result)
    double layer_factor = 1.4; // layer growth factor

    //start with 7 layers
    int nr_layers = 7;
    double coeff_sum = 0;
    for (int i = 0; i < nr_layers; i++)
        coeff_sum = coeff_sum + std::pow(layer_factor, i);

    double first_layer = skin_max * skin_factor / coeff_sum;

    // coefficient for 4 layers
    double coeff_4 = 0;
    for (int i = 0; i < 4; i++)
        coeff_4 = coeff_4 + std::pow(layer_factor, i);

    // increase number of layers while minimal thickness has less than 4 layers
    while (skin_min < first_layer * coeff_4)
    {
        nr_layers++;

        //recalculate first layer thickness
        coeff_sum = 0;
        for (int i = 0; i < nr_layers; i++)
            coeff_sum = coeff_sum + std::pow(layer_factor, i);

        first_layer = skin_max * skin_factor / coeff_sum;

        // Should this limitation be here?
        if (nr_layers >= 20)
            break;
    }

    SalomeVLParams svl;
    svl.total_thickness = 0;
    svl.nr_of_layers = nr_layers;
    svl.expansion_factor = layer_factor;
    for (int i = 0; i < nr_layers; i++)
        svl.total_thickness += first_layer * std::pow(layer_factor, i);

    return svl;
}

ShapeType MeshParameters::getType()
{
    return param_type;
}

std::vector<double> MeshParameters::getLayers()
{
    return thicknesses;
}