#include "geometryData.h"
#include <BRep_Builder.hxx>
#include <Bnd_Box.hxx>
#include <BRepBndLib.hxx>
#include <BRepTools.hxx>
#include <TopExp.hxx>
#include <TopExp_Explorer.hxx>
#include <future>
#include <thread>
#include "easylogging++/easylogging++.h"
#include "misc/miscFunctions.hpp"
#include "shapePartition.h"
#include "shapeAlgorithms.h"
#include "ShapeRepair.h"
#include "surroundingBuilder.h"
#include "cenos_exception.h"

#include <omp.h>

GeometryData::GeometryData(json input_json, std::string w)
{
    axisymmetric = false;
    changed_topology = false;
    wdir = w;

    // read operations than have to be performed before or during construction of geometryData
    readOperations(input_json);
    // if no geometry reading functions are passed, read json input
    if (read_functions.empty())
        readInput(input_json);
    else
    {
        for (auto fn : read_functions)
            fn();
    }

    // do operations that need to be performed on current geometry
    for (auto fn : operate_functions)
        fn();

    update();
}


GeometryData::GeometryData(std::string w)
{
    unit = LengthUnit::METER;
    wdir = w;
    dimension = 0;
    changed_topology = false;
    axisymmetric = false;
    editor_file = "";
}

GeometryData::GeometryData()
{
    wdir = "";
    unit = LengthUnit::METER;
    dimension = 0;
    changed_topology = false;
    axisymmetric = false;
    editor_file = "";
}

GeometryData::GeometryData(GeometryData& that)
{
    unit = that.unit;
    wdir = that.wdir;
    dimension = that.dimension;
    entities = that.entities;
    boundaries = that.boundaries;
    domains = that.domains;
    total_shape = that.total_shape;
    changed_topology = that.changed_topology;
    axisymmetric = that.axisymmetric;
    editor_file = that.editor_file;
}


void GeometryData::readInput(json input_json)
{
    if (input_json.is_null())
    {
        throw cenos_exception("Empty json object passed to GeometryData.");
    }

    std::vector<std::string> required_keys = { "solids", "faces", "edges" };

    for (auto& key : required_keys)
    {
        if (input_json.find(key) == input_json.end())
        {
            std::string message = "Error while reading geometry data: Could not find key " + key + " in geometryData json";
            throw cenos_exception(message);
        }
    }

    if (input_json["hasMeshData"].get<bool>())
    {
        meshfile = input_json["meshFile"].get<std::string>();
    }
    else
    {
        meshfile = "";
    }


    if (input_json.find("groupMode") != input_json.end())
        if (input_json.at("groupMode").get<std::string>() == "groupsFromPreprocessor")
            group_mode = GROUPS_FROM_PREPROCESSOR;
        else
            group_mode = GROUPS_FROM_ALL_ENTITIES;
    else
        group_mode = GROUPS_FROM_ALL_ENTITIES;


    if (input_json.find("geometryUnit") != input_json.end())
        unit = stringToLengthUnit(input_json.at("geometryUnit").get<std::string>());
    else
        unit = LengthUnit::NOT_SET;

    //read entities (solids, faces, edges)
    for (const auto& category : input_json.items())
    {
        if (category.key() == "solids" or
            category.key() == "faces" or
            category.key() == "edges")
        {
            // TO DO make reading entities in parallel
            for (const auto& item : category.value().items())
            {
                std::shared_ptr<TopoEntity> te = std::make_shared<TopoEntity>(item.value(), wdir);
                if (te->getId() == 0)
                    te->setId(newUniqueEntityId());
                entities.push_back(te);
            }
        }

        if (category.key() == "relations")
        {
            for (const auto& item : category.value().items())
            {
                relations[item.key()] = item.value().get<std::vector<std::string>>();
            }
        }
    }

    if (input_json.find("dimensions") != input_json.end())
        dimension = input_json.at("dimensions").get<int>();
    else
    {
        updateDimension();
    }

    for (const auto& category : input_json.items())
    {
        if (category.key() == "domains")
        {

            for (const auto& item : category.value().items())
            {
                std::shared_ptr<DomainGroup> gd = std::make_shared<DomainGroup>(item.value(), entities);
                domains.push_back(gd);
            }
        }

        if (category.key() == "boundaries")
        {

            for (const auto& item : category.value().items())
            {
                std::shared_ptr<BoundaryGroup> gb = std::make_shared<BoundaryGroup>(item.value(), entities);
                boundaries.push_back(gb);
            }
        }
    }

    if (input_json.find("partitionShape") != input_json.end())
    {
        if (!input_json.at("partitionShape").empty()) {
            std::shared_ptr<TopoEntity> te = std::make_shared<TopoEntity>(input_json.at("partitionShape"), wdir);
            total_shape = te->getShape();
        }
    }

    if (input_json.find("editorFile") != input_json.end())
    {
        editor_file = input_json.at("editorFile").get<std::string>();
    }
    else
        editor_file = "";

    // Reading input data done
}

void GeometryData::readOperations(json input_json)
{
    if (input_json.find("operate") != input_json.end())
    {

        for (const auto& op : input_json.at("operate").items())
        {
            json operate = op.value();

            if (operate.at("function") == "partition")
            {
                //skip partitioning if addAir call is found
                if (input_json.at("operate").find("partition") == input_json.at("operate").end())
                {
                    continue;
                }
                std::function<void()> partition_bind = std::bind(&GeometryData::partition, this);
                operate_functions.push_back(partition_bind);
            }
            else if (operate.at("function") == "addAir")
            {
                //pass zero padding, if it is not received as argument
                double padding = 0;
                std::vector<std::string> cut_boundaries;
                if (operate.find("padding") != operate.end())
                    padding = operate.at("padding").get<double>();
                if (operate.find("terminals") != operate.end())
                {
                    for (const auto& item : operate.at("terminals"))
                        cut_boundaries.push_back(item.at("name").get<std::string>());
                }
                std::function<void()> add_air_bind = std::bind(&GeometryData::addAir, this, padding, cut_boundaries);
                operate_functions.push_back(add_air_bind);
            }
            else if (operate.at("function") == "readCADfiles")
            {
                std::vector<std::string> filenames = operate.at("fileNames").get<std::vector<std::string>>();
                std::function<void()> readcads_bind = std::bind(&GeometryData::readCADfiles, this, filenames);
                read_functions.push_back(readcads_bind);
                // this group mode has to be used for FromCAD
                group_mode = GROUPS_FROM_ALL_ENTITIES;
            }
            else if (operate.at("function") == "scale")
            {
                double scale = operate.at("scale").get<double>();
                std::function<void()> scaleall_bind = std::bind(&GeometryData::scaleAll, this, scale);
                operate_functions.push_back(scaleall_bind);
            }
        }
    }

    if (input_json.find("variables") != input_json.end()) {
        variables = input_json["variables"];
    }


    //Operations done
}

void GeometryData::update()
{

    validate();
    // if geometry data is constructed manually, dimension is still zero!
    if (dimension == 0)
        updateDimension();

    if (dimension == 2)
        orientFaces();


    // Update domain and boundary groups that are passed in
    // Update necessary only if topology changing operations are performed
    if (changed_topology)
    {
        #pragma omp parallel  for
        for (int i = 0; i < domains.size(); i++)
        {
            domains[i]->clearEntities();
            domains[i]->clearBoundaries();
            for (auto ent : entities)
            {
                if (ent->getDimension() != dimension)
                    continue;
                if (ShapeAlgorithms::have_equal_shapes(ent->getShape(), domains[i]->getShape()))
                    domains[i]->addEntity(ent);
            }

            //reset id as after topology change indefinite boundaries/domains may have overlapping ids
            domains[i]->setId(0);
            domains[i]->update();
        }

        #pragma omp parallel  for
        for (int i = 0; i < boundaries.size(); i++)
        {
            boundaries[i]->clearEntities();
            boundaries[i]->clearDomains();
            for (auto ent : entities)
            {
                if (ent->getDimension() != dimension - 1)
                    continue;

                if (ShapeAlgorithms::have_equal_shapes(ent->getShape(), boundaries[i]->getShape()))
                    boundaries[i]->addEntity(ent);
            }
            //reset id as after topology change indefinite boundaries/domains may have overlapping ids
            boundaries[i]->setId(0);
            boundaries[i]->update();
        }

        //reset back, since next time update does not require domain and boundary group re-check
        changed_topology = false;
    }


    if (total_shape.IsNull())
    {
        ShapePartition part(entities);
        total_shape = part.getPartition();
    }

    //create indefinite boundary groups
    std::vector<std::string> used_face_names;
    for (auto bnd : this->getBoundaries())
    {
        for (auto ent : bnd->getEntities())
        {
            if (std::find(used_face_names.begin(), used_face_names.end(), ent->getName()) == used_face_names.end())
                used_face_names.push_back(ent->getName());
        }
    }

    for (auto ent : entities)
    {
        if (ent->getDimension() == dimension)
            continue;

        if (ent->getDimension() == dimension-2)
            continue;

        if (std::find(used_face_names.begin(), used_face_names.end(), ent->getName()) == used_face_names.end())
        {
            std::shared_ptr<BoundaryGroup> bgroup = std::make_shared<BoundaryGroup>();
            bgroup->setName(ent->getName());
            bgroup->setLabel(ent->getName());
            bgroup->setId(newUniqueGroupId());
            if (group_mode == GROUPS_FROM_PREPROCESSOR)
            {
                bgroup->setRole("indefinite");
                bgroup->setIndefinite(true);
            }
            else
            {
                bgroup->setRole("");
                bgroup->setIndefinite(false);
            }
            bgroup->addEntity(ent);
            bgroup->update();
            boundaries.push_back(bgroup);
        }
    }

    //create indefinite domain groups
    std::vector<std::string> used_solid_names;
    for (auto dom : this->getDomains())
    {
        for (auto ent : dom->getEntities())
        {
            if (std::find(used_solid_names.begin(), used_solid_names.end(), ent->getName()) == used_solid_names.end())
                used_solid_names.push_back(ent->getName());
        }
    }

    for (auto ent : entities)
    {
        if (ent->getDimension() != dimension)
            continue;

        if (std::find(used_solid_names.begin(), used_solid_names.end(), ent->getName()) == used_solid_names.end())
        {
            std::shared_ptr<DomainGroup> dgroup = std::make_shared<DomainGroup>();
            dgroup->setName(ent->getName());
            dgroup->setLabel(ent->getName());
            dgroup->setRole("");
            dgroup->setId(newUniqueGroupId());
            dgroup->addEntity(ent);
            dgroup->update();
            domains.push_back(dgroup);
        }
    }

    updateRelations();
}


void GeometryData::writeBreps() const
{
    misc::createFolder(wdir + "/geometry");
    misc::createFolder(wdir + "/geometry/raw/");

    #pragma omp parallel  for
    for (int i = 0; i < entities.size(); i++)
    {
        entities[i]->writeBrep(wdir);
    }
}


json GeometryData::getJson() const
{
    json out_json;
    out_json["solids"] = json::array();
    out_json["faces"] = json::array();
    out_json["edges"] = json::array();
    out_json["editorFile"] = editor_file;

    for (auto ent : entities)
    {
        if (ent->getType() == ShapeType::EDGE)
            out_json["edges"].push_back(ent->getJson());
        else if (ent->getType() == ShapeType::FACE)
            out_json["faces"].push_back(ent->getJson());
        else if (ent->getType() == ShapeType::SOLID)
            out_json["solids"].push_back(ent->getJson());
    }

    out_json["relations"] = relations;

    json domains_json;

    std::string vtk_preview_path = wdir + "/vtk";
    misc::createFolder(vtk_preview_path);

#pragma omp parallel for
    for (int i = 0; i < domains.size(); i++)
    {
        json s = domains[i]->getJson();

        MeshEntityStats stat;
        for (auto ent : domains[i]->getEntities())
        {
            if (mesh_stats.hasEntity(ent->getName()))
                stat = stat + mesh_stats.entity_stats.at(ent->getName());
        }
        s["meshStatistics"] = stat.getJson();

        //this is most time consuming part
        domains[i]->writePreviewVTK(wdir);

#pragma omp critical
        domains_json.push_back(s);
    }

    //loop over boundaries
    json boundaries_json = json::array();

    #pragma omp parallel for
    for (int i = 0; i< boundaries.size(); i++)
    {
        json b = boundaries[i]->getJson();

        MeshEntityStats stat;
        for (auto ent : boundaries[i]->getEntities())
        {
            if (mesh_stats.hasEntity(ent->getName()))
                stat = stat + mesh_stats.entity_stats.at(ent->getName());
        }
        b["meshStatistics"] = stat.getJson();

        //this is most time consuming part
        boundaries[i]->writePreviewVTK(wdir);

#pragma omp critical
        boundaries_json.push_back(b);
    }

    out_json["domains"] = domains_json;
    out_json["boundaries"] = boundaries_json;

    ShapeStats stat;
    //initialize stat with values from first entity. Otherwise default values (zeros) may remain through whole entities
    ShapeAlgorithms::getShapeStatistics(&entities[0]->getShape(), stat);
    for (auto ent : entities)
    {
        if (ent->getDimension() != dimension)
            continue;

        ShapeStats entity_stat;
        ShapeAlgorithms::getShapeStatistics(&ent->getShape(), entity_stat);
        stat = stat + entity_stat;
    }

    out_json["statistics"] = stat.getJson();

    out_json["meshStatistics"] = mesh_stats.total_stats.getJson();

    if (group_mode == GROUPS_FROM_ALL_ENTITIES)
        out_json["groupMode"] = "groupsFromAllEntities";
    else if (group_mode == GROUPS_FROM_PREPROCESSOR)
        out_json["groupMode"] = "groupsFromPreprocessor";

    out_json["dimensions"] = dimension;
    out_json["geometryUnit"] = lengthUnitToString(unit);

    std::shared_ptr<TopoEntity> te = std::make_shared<TopoEntity>(total_shape, "partitionShape");
    te->writeBrep(wdir);
    out_json["partitionShape"] = te->getJson();

    if (meshfile.length() > 0)
    {
        out_json["meshFile"] = meshfile;
        out_json["hasMeshData"] = true;
    }
    else
    {
        checkIfMeshable();
        out_json["meshFile"] = json();
        out_json["hasMeshData"] = false;
    }

    out_json["variables"] = variables;

    return out_json;
}


LengthUnit GeometryData::getLengthUnit()
{
    return unit;
}


void GeometryData::setMeshFile(std::string mf)
{
    meshfile = mf;
}

void GeometryData::setMeshStats(MeshStats ms)
{
    mesh_stats = ms;
}


std::vector<std::shared_ptr<DomainGroup>> GeometryData::getDomains() const
{
    return domains;
}


std::vector<std::shared_ptr<BoundaryGroup>> GeometryData::getBoundaries() const
{
    return boundaries;
}

std::vector<std::shared_ptr<TopoEntity>> GeometryData::getEntities() const
{
    return entities;
}

std::shared_ptr<DomainGroup> GeometryData::getDomainByName(std::string n)
{
    for (auto dom : domains)
    {
        if (dom->getName() == n)
            return dom;
    }
    std::string message = "Error in geometry data: Could not find domain " + n + ". Check your geometry input!";
    throw cenos_exception(message);
}

std::shared_ptr<BoundaryGroup> GeometryData::getBoundaryByName(std::string n)
{
    for (auto bnd : boundaries)
    {
        if (bnd->getName() == n)
            return bnd;
    }
    std::string message = "Error in geometry data: Could not find boundary " + n + ". Check your geometry input!";
    throw cenos_exception(message);
}

std::shared_ptr<TopoEntity> GeometryData::getEntityByName(std::string n)
{
    for (auto ent : entities)
    {
        if (ent->getName() == n)
            return ent;
    }
    std::string message = "Error in geometry data: Could not find entity " + n + ". Check your geometry input!";
    throw cenos_exception(message);
}

void GeometryData::addEntity(std::shared_ptr<TopoEntity> e)
{
    if (e->getId() == 0)
        e->setId(newUniqueEntityId());
    entities.push_back(e);
}

void GeometryData::addDomainGroup(std::shared_ptr<DomainGroup> d)
{
    domains.push_back(d);
}

void GeometryData::addBoundaryGroup(std::shared_ptr<BoundaryGroup> b)
{
    if (b->getId() == 0)
        b->setId(newUniqueGroupId());
    boundaries.push_back(b);
}



void GeometryData::partition()
{
    ShapePartition part(entities);

    TopoDS_Shape partition_shape = part.getPartition();
    this->entities.clear();

    TopTools_IndexedMapOfShape emap;
    TopExp::MapShapes(partition_shape, TopAbs_EDGE, emap);
    for (int ei = 1; ei <= emap.Extent(); ei++)
    {
        std::shared_ptr<TopoEntity> te = std::make_shared<TopoEntity>(emap.FindKey(ei), "Edge_" + std::to_string(ei), ei);
        this->addEntity(te);
    }

    TopTools_IndexedMapOfShape fmap;
    TopExp::MapShapes(partition_shape, TopAbs_FACE, fmap);
    for (int fi = 1; fi <= fmap.Extent(); fi++)
    {
        std::shared_ptr<TopoEntity> te = std::make_shared<TopoEntity>(fmap.FindKey(fi), "Face_" + std::to_string(fi), fi + emap.Extent());
        this->addEntity(te);
    }

    TopTools_IndexedMapOfShape smap;
    TopExp::MapShapes(partition_shape, TopAbs_SOLID, smap);
    for (int si = 1; si <= smap.Extent(); si++)
    {
        std::shared_ptr<TopoEntity> te = std::make_shared<TopoEntity>(smap.FindKey(si), "Solid_" + std::to_string(si), si + emap.Extent() + fmap.Extent());
        this->addEntity(te);
    }

    for (int fi = 1; fi <= fmap.Extent(); fi++)
    {
        relations["Face_" + std::to_string(fi)] = {};
    }

    for (int si = 1; si <= smap.Extent(); si++)
    {
        relations["Solid_" + std::to_string(si)] = {};

        for (TopExp_Explorer face_explorer(smap.FindKey(si), TopAbs_ShapeEnum::TopAbs_FACE); face_explorer.More(); face_explorer.Next())
        {
            for (int fi = 1; fi <= fmap.Extent(); fi++)
            {
                if (face_explorer.Current().IsSame(fmap.FindKey(fi)))
                {
                    relations["Solid_" + std::to_string(si)].push_back("Face_" + std::to_string(fi));
                    relations["Face_" + std::to_string(fi)].push_back("Solid_" + std::to_string(si));
                }
            }
        }
    }

    changed_topology = true;
    this->writeBreps();
}

void GeometryData::addAir(double padding, std::vector<std::string> cut_boundary_names )
{
    if (dimension == 0)
        updateDimension();

    TopTools_ListOfShape cut_boundaries; // boundaries used to cut air box
    for (auto ent : entities)
    {
        if (std::find(cut_boundary_names.begin(), cut_boundary_names.end(), ent->getName()) != cut_boundary_names.end())
            cut_boundaries.Append(ent->getShape());
    }

    SurroundingBuilder sbuilder;
    sbuilder.setCutBoundaries(cut_boundaries);
    sbuilder.setEntities(entities);
    sbuilder.setPadding(padding);

    if (dimension == 3)
        sbuilder.setType(SurroundingBuilder::FULL);
    else if (dimension == 2 and axisymmetric == false)
        sbuilder.setType(SurroundingBuilder::PLANAR_2D);
    else if (dimension == 2 and axisymmetric == true)
        sbuilder.setType(SurroundingBuilder::AXYSIMMETRIC_2D);
    else
        throw(cenos_exception("Dimension of GeometryData is " + std::to_string(dimension) + "is not supported"));

    sbuilder.build();  // this command may change entities!!!

    TopoDS_Shape air_infinity = sbuilder.getInfinity();
    TopoDS_Shape air_shape = sbuilder.getSurrounding();
    TopoDS_Shape air_symmetry = sbuilder.getSymmetry(); // currently only for 2D axisymmetric

    ShapePartition new_part(entities);

    if (dimension == 3)
    {
        for (TopExp_Explorer shape_explorer(air_shape, TopAbs_SOLID); shape_explorer.More(); shape_explorer.Next())
            new_part.addShape(shape_explorer.Current());

    }
    if (dimension == 2)
    {
        for (TopExp_Explorer shape_explorer(air_shape, TopAbs_FACE); shape_explorer.More(); shape_explorer.Next())
            new_part.addShape(shape_explorer.Current());
    }

    TopoDS_Shape new_partition = new_part.getPartition();
    total_shape = new_partition;

    this->entities.clear();
    if (group_mode == Mode::GROUPS_FROM_ALL_ENTITIES)
    {
        this->domains.clear();
        this->boundaries.clear();
    }
    this->relations.clear();

    TopTools_IndexedMapOfShape emap;
    TopExp::MapShapes(new_partition, TopAbs_EDGE, emap);
    for (int ei = 1; ei <= emap.Extent(); ei++)
    {
        std::shared_ptr<TopoEntity> te = std::make_shared<TopoEntity>(emap.FindKey(ei), "Edge_" + std::to_string(ei), ei);
        this->addEntity(te);
    }

    TopTools_IndexedMapOfShape fmap;
    TopExp::MapShapes(new_partition, TopAbs_FACE, fmap);
    for (int fi = 1; fi <= fmap.Extent(); fi++)
    {
        std::shared_ptr<TopoEntity> te = std::make_shared<TopoEntity>(fmap.FindKey(fi), "Face_" + std::to_string(fi ), fi + emap.Extent());
        this->addEntity(te);
    }

    TopTools_IndexedMapOfShape smap;
    TopExp::MapShapes(new_partition, TopAbs_SOLID, smap);
    for (int si = 1; si <= smap.Extent(); si++)
    {
        std::shared_ptr<TopoEntity> te = std::make_shared<TopoEntity>(smap.FindKey(si), "Solid_" + std::to_string(si ), si + emap.Extent() + fmap.Extent() );
        this->addEntity(te);

    }

    std::shared_ptr<DomainGroup> dgroup = std::make_shared<DomainGroup>();
    dgroup->setName("air");
    dgroup->setLabel("Air");
    dgroup->setRole("air");
    for (auto ent : entities)
    {
        if (ent->getDimension() != dimension)
            continue;

        if (ShapeAlgorithms::have_equal_shapes(ent->getShape(), air_shape))
        {
            dgroup->setId(ent->getId());
            dgroup->addEntity(ent);
        }
    }
    dgroup->update();
    domains.push_back(dgroup);

    if (!air_infinity.IsNull())
    {
        std::shared_ptr<BoundaryGroup> bgroup = std::make_shared<BoundaryGroup>();
        bgroup->setName("air_infinity");
        bgroup->setLabel("Infinity");
        bgroup->setRole("infinity");
        for (auto ent : entities)
        {
            if (ent->getDimension() != dimension - 1)
                continue;

            if (ShapeAlgorithms::have_equal_shapes(ent->getShape(), air_infinity))
            {
                bgroup->setId(ent->getId());
                bgroup->addEntity(ent);
            }
        }
        bgroup->setIndefinite(false);
        bgroup->update();
        boundaries.push_back(bgroup);
    }

    if (misc::getCenosApp() == "INDUCTION" and dimension == 3) {
        // Lets add the terminals as boundaries and set their role, 
        // so the user doesnt have to do it himself.

        int terminal_number = 1;
        for (TopoDS_Shape terminal : cut_boundaries) {
            if (!terminal.IsNull())
            {
                std::shared_ptr<BoundaryGroup> bgroup = std::make_shared<BoundaryGroup>();
                for (auto ent : entities)
                {
                    if (ent->getDimension() != dimension - 1)
                        continue;

                    if (ShapeAlgorithms::have_equal_shapes(ent->getShape(), terminal))
                    {
                        bgroup->setId(ent->getId());
                        bgroup->addEntity(ent);
                    }
                }
                bgroup->setIndefinite(false);
                bgroup->update();
                boundaries.push_back(bgroup);

                std::string terminalName = "Terminal";
                std::string terminalLabel = "Terminal";
                std::string terminal_index = std::to_string(terminal_number);


                bgroup->setName(terminalName + terminal_index);
                bgroup->setLabel(terminalLabel + terminal_index);
                if (terminal_number % 2 != 0) {
                    bgroup->setRole("terminal");
                }
                else {
                    bgroup->setRole("terminal2");
                }
                terminal_number++;
            }
        }
    }

    if (!air_symmetry.IsNull())
    {
        std::shared_ptr<BoundaryGroup> bgroup = std::make_shared<BoundaryGroup>();
        bgroup->setName("air_axis");
        bgroup->setRole("axis");
        bgroup->setLabel("Symmetry Axis");
        for (auto ent : entities)
        {
            if (ent->getDimension() != dimension - 1)
                continue;

            if (ShapeAlgorithms::have_equal_shapes(ent->getShape(), air_symmetry))
            {
                bgroup->setId(ent->getId());
                bgroup->addEntity(ent);
            }
        }
        bgroup->setIndefinite(false);
        bgroup->update();
        boundaries.push_back(bgroup);
    }


    if (dimension == 2)
    {
        for (int ei = 1; ei <= emap.Extent(); ei++)
        {
            relations["Edge_" + std::to_string(ei)] = {};
        }

        for (int fi = 1; fi <= fmap.Extent(); fi++)
        {
            relations["Face_" + std::to_string(fi)] = {};

            for (TopExp_Explorer edge_explorer(fmap.FindKey(fi), TopAbs_ShapeEnum::TopAbs_EDGE); edge_explorer.More(); edge_explorer.Next())
            {
                for (int ei = 1; ei <= emap.Extent(); ei++)
                {
                    if (edge_explorer.Current().IsSame(emap.FindKey(ei)))
                    {
                        relations["Face_" + std::to_string(fi)].push_back("Edge_" + std::to_string(ei));
                        relations["Edge_" + std::to_string(ei)].push_back("Face_" + std::to_string(fi));
                    }
                }
            }
        }
    }
    else if (dimension == 3)
    {
        for (int fi = 1; fi <= fmap.Extent(); fi++)
        {
            relations["Face_" + std::to_string(fi)] = {};
        }

        for (int si = 1; si <= smap.Extent(); si++)
        {
            relations["Solid_" + std::to_string(si)] = {};

            for (TopExp_Explorer face_explorer(smap.FindKey(si), TopAbs_ShapeEnum::TopAbs_FACE); face_explorer.More(); face_explorer.Next())
            {
                for (int fi = 1; fi <= fmap.Extent(); fi++)
                {
                    if (face_explorer.Current().IsSame(fmap.FindKey(fi)))
                    {
                        relations["Solid_" + std::to_string(si)].push_back("Face_" + std::to_string(fi));
                        relations["Face_" + std::to_string(fi)].push_back("Solid_" + std::to_string(si));
                    }
                }
            }
        }
    }


    changed_topology = true;
    this->writeBreps();
}

void GeometryData::merge(std::vector<std::string> group_names, std::string label)
{
    std::vector<std::shared_ptr<Group>> to_merge_groups;
    bool is_bounary = false;
    for (auto bnd : boundaries)
    {
        bnd->setPreviewWriting(false);
        if (std::find(group_names.begin(), group_names.end(), bnd->getName()) != group_names.end())
        {
            to_merge_groups.push_back(bnd);
            is_bounary = true;
        }
    }

    for (auto dom : domains)
    {
        dom->setPreviewWriting(false);
        if (std::find(group_names.begin(), group_names.end(), dom->getName()) != group_names.end())
            to_merge_groups.push_back(dom);
    }

    json gr_json;
    gr_json["name"] = "Merged_" + to_merge_groups[0]->getName();
    gr_json["label"] = label;
    for (auto gr : to_merge_groups)
    {
        // what to do with roles of merged groups?
        // bgr->setRole("");

        gr_json["mergedFrom"].push_back(gr->getJson());
        for (auto ent : gr->getEntities())
            gr_json["entities"].push_back(ent->getName());
    }
    gr_json["role"] = "";
    gr_json["id"] = newUniqueGroupId();

    if (is_bounary)
    {
        //insert merged boundaries
        auto it = boundaries.begin();
        while (it != boundaries.end())
        {
            if (std::find(to_merge_groups.begin(), to_merge_groups.end(), *it) != to_merge_groups.end())
            {
                std::shared_ptr<BoundaryGroup> bgr = std::make_shared<BoundaryGroup>(gr_json, entities);
                boundaries.insert(it, bgr);
                break;
            }
            else
                ++it;
        }

        //remove old boundaries
        it = boundaries.begin();
        while (it != boundaries.end())
        {
            if (std::find(to_merge_groups.begin(), to_merge_groups.end(), *it) != to_merge_groups.end())
                it = boundaries.erase(it);
            else
                ++it;
        }

    }
    else
    {
        //insert merged domains
        auto it = domains.begin();
        while (it != domains.end())
        {
            if (std::find(to_merge_groups.begin(), to_merge_groups.end(), *it) != to_merge_groups.end())
            {
                std::shared_ptr<DomainGroup> dgr = std::make_shared<DomainGroup>(gr_json, entities);
                domains.insert(it, dgr);
                break;
            }
            else
                ++it;
        }

        //remove old domains
        it = domains.begin();
        while (it != domains.end())
        {
            if (std::find(to_merge_groups.begin(), to_merge_groups.end(), *it) != to_merge_groups.end())
                it = domains.erase(it);
            else
                ++it;
        }
    }

    for (auto bnd : boundaries)
        bnd->clearDomains();

    for (auto dom : domains)
        dom->clearBoundaries();

    updateRelations();
}



void GeometryData::unmerge(std::vector<std::string> group_names)
{
    for (auto group_name : group_names)
    {

        bool is_bounary = false;
        std::shared_ptr<Group> group;
        for (auto bnd : boundaries)
        {
            bnd->setPreviewWriting(false);
            if (group_name== bnd->getName())
            {
                group = bnd;
                is_bounary = true;
                break;
            }

        }

        for (auto dom : domains)
        {
            dom->setPreviewWriting(false);
            if (group_name== dom->getName())
            {
                group = dom;
                break;
            }
        }


        if (is_bounary)
        {
            //insert merged boundaries
            auto it = boundaries.begin();
            while (it != boundaries.end())
            {
                if (group == *it)
                {
                    std::vector<std::shared_ptr<BoundaryGroup>> unmerged_boundaries;
                    for (auto gr : group->getMergedFrom())
                    {
                        std::shared_ptr<BoundaryGroup> bgr = std::make_shared<BoundaryGroup>(gr->getJson(), entities);
                        unmerged_boundaries.push_back(bgr);
                    }
                    boundaries.insert(it, unmerged_boundaries.begin(), unmerged_boundaries.end());
                    break;
                }
                else
                    ++it;
            }

            //remove old boundaries
            it = boundaries.begin();
            while (it != boundaries.end())
            {
                if (group == *it)
                    it = boundaries.erase(it);
                else
                    ++it;
            }

        }
        else
        {
            //insert merged domains
            auto it = domains.begin();
            while (it != domains.end())
            {
                if (group == *it)
                {
                    std::vector< std::shared_ptr<DomainGroup> > unmerged_domains;

                    for (auto gr : group->getMergedFrom())
                    {
                        std::shared_ptr<DomainGroup> dgr = std::make_shared<DomainGroup>(gr->getJson(), entities);
                        unmerged_domains.push_back(dgr);
                    }

                    domains.insert(it, unmerged_domains.begin(), unmerged_domains.end());
                    break;
                }
                else
                    ++it;
            }

            //remove old domains
            it = domains.begin();
            while (it != domains.end())
            {
                if (group == *it)
                    it = domains.erase(it);
                else
                    ++it;
            }
        }
    }

    for (auto bnd : boundaries)
        bnd->clearDomains();

    for (auto dom : domains)
        dom->clearBoundaries();

    updateRelations();
}


// Fixme: Use parameter as reference
void GeometryData::readCADfiles(std::vector<std::string> filenames)
{
    ShapePartition part;
    for (auto file_name : filenames)
    {
        std::string full_path = wdir + "/" + file_name;
        std::string ext = misc::getPathExtension(full_path);
        // make extension upper case. cad files can have mixed case symbols in extension
        std::transform(ext.begin(), ext.end(), ext.begin(), ::toupper);

        TopoDS_Shape shape;
        if (ext == ".STP" or ext == ".STEP")
        {
            shape = ShapeAlgorithms::readStepFile(full_path);
        }
        else if (ext == ".IGS" or ext == ".IGES")
        {
            shape = ShapeAlgorithms::readIgesFile(full_path);
        }
        else
        {
            throw(cenos_exception("File format " + ext + " is not supported. Please choose other CAD file!"));
        }
        shape = ShapeRepair::pre_process_cad(shape);
         
        TopAbs_ShapeEnum shapetype = ShapeAlgorithms::getHighestShapeType(shape);

        TopTools_IndexedMapOfShape shape_map;
        TopExp::MapShapes(shape, shapetype, shape_map);
        for (int i = 1; i <= shape_map.Extent(); i++)
            part.addShape(shape_map.FindKey(i));
    }

    TopoDS_Shape partition_shape = part.getPartition();
    total_shape = partition_shape;

    TopTools_IndexedMapOfShape emap;
    TopExp::MapShapes(partition_shape, TopAbs_EDGE, emap);
    for (int ei = 1; ei <= emap.Extent(); ei++)
    {
        std::shared_ptr<TopoEntity> te = std::make_shared<TopoEntity>(emap.FindKey(ei), "Edge_" + std::to_string(ei), ei);
        this->addEntity(te);
        dimension = 1;
    }

    TopTools_IndexedMapOfShape fmap;
    TopExp::MapShapes(partition_shape, TopAbs_FACE, fmap);
    for (int fi = 1; fi <= fmap.Extent(); fi++)
    {
        std::shared_ptr<TopoEntity> te = std::make_shared<TopoEntity>(fmap.FindKey(fi), "Face_" + std::to_string(fi), fi + emap.Extent());
        this->addEntity(te);
        dimension = 2;
    }

    TopTools_IndexedMapOfShape smap;
    TopExp::MapShapes(partition_shape, TopAbs_SOLID, smap);
    for (int si = 1; si <= smap.Extent(); si++)
    {
        std::shared_ptr<TopoEntity> te = std::make_shared<TopoEntity>(smap.FindKey(si), "Solid_" + std::to_string(si), si + emap.Extent() + fmap.Extent());
        this->addEntity(te);
        dimension = 3;
    }

    if (dimension == 2)
    {
        for (int ei = 1; ei <= emap.Extent(); ei++)
        {
            relations["Edge_" + std::to_string(ei)] = {};
        }

        for (int fi = 1; fi <= fmap.Extent(); fi++)
        {
            relations["Face_" + std::to_string(fi)] = {};

            for (TopExp_Explorer edge_explorer(fmap.FindKey(fi), TopAbs_ShapeEnum::TopAbs_EDGE); edge_explorer.More(); edge_explorer.Next())
            {
                for (int ei = 1; ei <= emap.Extent(); ei++)
                {
                    if (edge_explorer.Current().IsSame(emap.FindKey(ei)))
                    {
                        relations["Face_" + std::to_string(fi)].push_back("Edge_" + std::to_string(ei));
                        relations["Edge_" + std::to_string(ei)].push_back("Face_" + std::to_string(fi));
                    }
                }
            }
        }
    }
    else if (dimension == 3)
    {
        for (int fi = 1; fi <= fmap.Extent(); fi++)
        {
            relations["Face_" + std::to_string(fi)] = {};
        }

        for (int si = 1; si <= smap.Extent(); si++)
        {
            relations["Solid_" + std::to_string(si)] = {};

            for (TopExp_Explorer face_explorer(smap.FindKey(si), TopAbs_ShapeEnum::TopAbs_FACE); face_explorer.More(); face_explorer.Next())
            {
                for (int fi = 1; fi <= fmap.Extent(); fi++)
                {
                    if (face_explorer.Current().IsSame(fmap.FindKey(fi)))
                    {
                        relations["Solid_" + std::to_string(si)].push_back("Face_" + std::to_string(fi));
                        relations["Face_" + std::to_string(fi)].push_back("Solid_" + std::to_string(si));
                    }
                }
            }
        }
    }

    changed_topology = true;
    this->writeBreps();
}

void GeometryData::scaleAll(double scale)
{
    for (auto ent : entities)
    {
        if (!ent->isEnabled())
            continue;
        TopoDS_Shape sh = ent->getShape();
        ent->setShape(ShapeAlgorithms::scaleShape(sh, scale));
    }

    writeBreps();

    for (auto bnd : boundaries)
        bnd->update();

    for (auto dom : domains)
        dom->update();

    recalculateShape();
}


int GeometryData::newUniqueGroupId()
{
    int max_id = 0;
    for (auto bnd : boundaries)
    {
        if (max_id < bnd->getId())
            max_id = bnd->getId();
    }

    for (auto dom : domains)
    {
        if (max_id < dom->getId())
            max_id = dom->getId();
    }
    return max_id + 1;
}

int GeometryData::newUniqueEntityId()
{
    int max_id = 0;
    for (auto ent : entities)
    {
        if (max_id < ent->getId())
            max_id = ent->getId();
    }

    return max_id + 1;
}
// updates relations object
// that keeps information about which face belongs to which solid
void GeometryData::updateRelations()
{
    if (relations.empty())
    {
        //initialize all vectors in entity map
        for (auto ent : entities)
        {
            relations[ent->getName()] = {};
        }

        //check all solids (faces in 2D) if they have face (edge in 2D) as subshape
        for (auto s_ent : entities)
        {
            if (s_ent->getDimension() != dimension)
                continue;

            for (auto f_ent : entities)
            {
                if (f_ent->getDimension() != dimension - 1)
                    continue;

                if (relations[f_ent->getName()].size() == 2)
                    continue;

                if (dimension == 3 and ShapeAlgorithms::have_common_faces(s_ent->getShape(), f_ent->getShape()))
                {
                    relations[s_ent->getName()].push_back(f_ent->getName());
                    relations[f_ent->getName()].push_back(s_ent->getName());
                }
                else if (dimension == 2 and ShapeAlgorithms::have_common_edges(s_ent->getShape(), f_ent->getShape()))
                {
                    relations[s_ent->getName()].push_back(f_ent->getName());
                    relations[f_ent->getName()].push_back(s_ent->getName());
                }
            }
        }
    }

    //We clear inDomain flag for every entity
    for (auto boundary : this->getBoundaries())
    {
        boundary->unsetIsindomain();

        for (auto bnd_ent : boundary->getEntities())
            bnd_ent->unsetIsindomain();
    }

    //std::cout << "INNER BOUNDARIES CLEARED" << std::endl;

    //update relations between domains and boundaries
    for (auto domain : this->getDomains())
    {
        //*************** CHANGE 1 *********
        //vector for domain entity names
        std::vector<std::string> domain_entity_names;
        for (auto dom_ent : domain->getEntities())
            domain_entity_names.push_back(dom_ent->getName());
        //***************
        for (auto dom_ent : domain->getEntities())
        {
            // loop over boundaries and their entities

            for (auto boundary : this->getBoundaries())
            {
                int flag_count = 0;

                
                for (auto bnd_ent : boundary->getEntities())
                {
                    //*************** CHANGE 2 *********
                    auto w = relations[bnd_ent->getName()];
                    // if w size is 1, it is external
                    // if size is 2, internal
                    if (w.size() == 2)
                    {
                        //now determine if this entity is between solids (faces in 2D) of same domain
                        //determine if both w[0] and w[1] belong to domain. Then it means that entity bnd_ent is between these two solids
                        if ((std::find(domain_entity_names.begin(), domain_entity_names.end(), w[0]) != domain_entity_names.end()) and
                            (std::find(domain_entity_names.begin(), domain_entity_names.end(), w[1]) != domain_entity_names.end()) and
                            !bnd_ent->getIsindomain())
                        {

                            bnd_ent->setIsindomain();
                            flag_count++;
                            //std::cout << "INNER BOUNDARY FLAG: [" << bnd_ent->getName() << "]" << std::endl;
                        }
                    }
                    //***************
                    auto v = relations[dom_ent->getName()];
                    if (std::find(v.begin(), v.end(), bnd_ent->getName()) != v.end())
                    {
                        domain->setBoundary(boundary);
                        boundary->setDomain(domain);
                    }
                }


                if (boundary->getEntities().size() == flag_count && flag_count > 0 && !boundary->getIsindomain())
                {
                    boundary->setIsindomain();
                   // std::cout << "BOUNDARY DOMAIN FLAGGED: [" << boundary->getName() << "]" << std::endl;
                }
             
            }




        }
    }

}

void GeometryData::updateDimension()
{
    // reset dimension to zero, then loop over all entities, assign highest of entity dimensions
    dimension = 0;
    for (auto ent : entities)
    {
        if (ent->getDimension() == 3)
        {
            dimension = 3;
            break;
        }
        else if (ent->getDimension() == 2 or ent->getDimension() == 1)
        {
            if (dimension < ent->getDimension())
                dimension = ent->getDimension();
        }
    }
}


// orient faces for 2D case that they all have same orientation
// if orientation changed, also recalculate partition
void GeometryData::orientFaces()
{
    bool orientation_changed = false;
    for (auto ent : entities)
    {
        if (ent->getDimension() == 2)
        {
            gp_Dir dir_ent;
            TopoDS_Shape sh = ent->getShape();
            dir_ent = ShapeAlgorithms::getFaceNormal(sh);

            // For some reason the surface normal for "Reversed" faces is not reversed.
            // So these "Reversed" faces are reversed to get "Forward" orientation.
            if (sh.Orientation() == TopAbs_Orientation::TopAbs_REVERSED) {
                sh.Reverse();
            }

            if (dir_ent.IsOpposite(gp_Dir(0, 0, 1), 0.05))
            {
                ent->setShape(sh.Reversed());
                orientation_changed = true;
            }
        }
    }
    if (orientation_changed)
        recalculateShape();
}


void GeometryData::setLengthUnit(LengthUnit u)
{
    unit = u;
}

void GeometryData::setGroupMode(Mode m)
{
    group_mode = m;
}

void GeometryData::validate() const
{
    // This error should never occur,
    // if this happens there is a critical error
    if (wdir == "")
        throw cenos_exception("Working directory not set in GeometryData. Please report this to Cenos support!");

    std::vector<std::string> group_names;

    for (auto dom : domains)
        group_names.push_back(dom->getName());

    for (auto bnd : boundaries)
        group_names.push_back(bnd->getName());

    std::sort(group_names.begin(), group_names.end());

    auto it = std::unique(group_names.begin(), group_names.end());
    // UI does the check, but for safety we check also here
    if (it != group_names.end())
        throw cenos_exception("Names of boundaries and domains are not unique! Please rename boundaries or domains!");
}



void GeometryData::checkIfMeshable() const
{
    //check iternal faces
    //currently netgen mesher cannot handle such faces
    for (auto s : entities)
    {
        if (s->getType() == ShapeType::SOLID)
        {
            if (ShapeAlgorithms::hasInternalFaces(s->getShape()))
            {
                std::string domname = "";
                for (auto dom : domains)
                {
                    auto ents = dom->getEntities();
                    if (std::find_if(ents.begin(), ents.end(), [s]( const std::shared_ptr<TopoEntity> gr) { return gr->getName()== s->getName() ; }) != ents.end())
                    {
                        domname = dom->getName();
                    }
                }

                std::string message = "WARNING: Currently Cenos MeshGenerator does not support internal faces. Internal face found in domain " + domname;
                std::cerr << message << std::endl;
            }
        }
    }
}

void GeometryData::setAxisymmetric(bool b)
{
    axisymmetric = b;
}

bool GeometryData::isAxisymmetric()
{
    return axisymmetric;
}

ShapeStats GeometryData::getShapeStatistics()
{
    ShapeStats stat;
    ShapeAlgorithms::getShapeStatistics(&getShape(), stat);
    return stat;
}


TopoDS_Shape GeometryData::getShape()
{
    if (total_shape.IsNull())
    {
        ShapePartition part(entities);
        total_shape = part.getPartition();
    }

    return total_shape;
}

void GeometryData::recalculateShape()
{
    // partition ( BRepAlgoAPI_BuilderAlgo ) produces bad results sometimes for 3D geometries that do not overlap (only touch)
    // therefore  GEOMAlgo_Gluer2  used for 3D, and BRepAlgoAPI_BuilderAlgo for 2D cases
    ShapePartition part(entities);
    if (dimension == 3)
        total_shape = part.getGlue();
    else if (dimension == 2)
        total_shape = part.getPartition();
    else
        throw(cenos_exception("Cannot calculate partition for geometryData with dimension other than 2 or 3."));
}

std::map<std::string, std::vector<std::string>> GeometryData::getRelations()
{
    return relations;
}

int GeometryData::getDimension()
{
    if (dimension == 0)
        updateDimension();
    return dimension;
}