ControlDeformExtSurface.cpp

/*---------------------------------------------------------------------------*\
 *
 *  mimmo
 *
 *  Copyright (C) 2015-2021 OPTIMAD engineering Srl
 *
 *  -------------------------------------------------------------------------
 *  License
 *  This file is part of mimmo.
 *
 *  mimmo is free software: you can redistribute it and/or modify it
 *  under the terms of the GNU Lesser General Public License v3 (LGPL)
 *  as published by the Free Software Foundation.
 *
 *  mimmo 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 Lesser General Public
 *  License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with mimmo. If not, see <http://www.gnu.org/licenses/>.
 *
\*---------------------------------------------------------------------------*/
#include "ControlDeformExtSurface.hpp"
#include "SkdTreeUtils.hpp"
#include <volcartesian.hpp>
#include <CG.hpp>
 
namespace mimmo{
 
ControlDeformExtSurface::ControlDeformExtSurface(){
    m_name = "mimmo.ControlDeformExtSurface";
    m_tolerance = 0.0;
 
    m_allowed.insert((FileType::_from_string("STL"))._to_integral());
    m_allowed.insert((FileType::_from_string("SURFVTU"))._to_integral());
    m_allowed.insert((FileType::_from_string("NAS"))._to_integral());
 
};
 
ControlDeformExtSurface::ControlDeformExtSurface(const bitpit::Config::Section & rootXML){
 
    m_name = "mimmo.ControlDeformExtSurface";
    m_tolerance = 0.0;
    m_allowed.insert((FileType::_from_string("STL"))._to_integral());
    m_allowed.insert((FileType::_from_string("SURFVTU"))._to_integral());
    m_allowed.insert((FileType::_from_string("NAS"))._to_integral());
 
    std::string fallback_name = "ClassNONE";
    std::string input = rootXML.get("ClassName", fallback_name);
    input = bitpit::utils::string::trim(input);
    if(input == "mimmo.ControlDeformExtSurface"){
        absorbSectionXML(rootXML);
    }else{
        warningXML(m_log, m_name);
    };
}
 
ControlDeformExtSurface::~ControlDeformExtSurface(){};
 
ControlDeformExtSurface::ControlDeformExtSurface(const ControlDeformExtSurface & other):BaseManipulation(other){
    m_allowed = other.m_allowed;
    m_geoList = other.m_geoList;
    m_geoFileList = other.m_geoFileList;
    m_tolerance = other.m_tolerance;
};
 
ControlDeformExtSurface & ControlDeformExtSurface::operator=(ControlDeformExtSurface other){
    swap(other);
    return(*this);
};
 
void ControlDeformExtSurface::swap(ControlDeformExtSurface & x) noexcept
{
    std::swap(m_allowed, x.m_allowed);
    std::swap(m_geoList, x.m_geoList);
    std::swap(m_geoFileList, x.m_geoFileList);
    std::swap(m_tolerance, x.m_tolerance);
    m_violationField.swap(x.m_violationField);
    m_defField.swap(x.m_defField);
    BaseManipulation::swap(x);
}
void
ControlDeformExtSurface::buildPorts(){
    bool built = true;
 
    built = (built && createPortIn<dmpvecarr3E*, ControlDeformExtSurface>(this, &mimmo::ControlDeformExtSurface::setDefField, M_GDISPLS));
    built = (built && createPortIn<MimmoSharedPointer<MimmoObject>, ControlDeformExtSurface>(this, &mimmo::ControlDeformExtSurface::setGeometry, M_GEOM, true));
    built = (built && createPortIn<MimmoSharedPointer<MimmoObject>, ControlDeformExtSurface>(this, &mimmo::ControlDeformExtSurface::addConstraint, M_GEOM2));
 
    built = (built && createPortOut<double, ControlDeformExtSurface>(this, &mimmo::ControlDeformExtSurface::getViolation, M_VALUED));
    built = (built && createPortOut<dmpvector1D*, ControlDeformExtSurface>(this, &mimmo::ControlDeformExtSurface::getViolationField, M_SCALARFIELD));
    m_arePortsBuilt = built;
};
 
double
ControlDeformExtSurface::getViolation(){
 
    double result = -1.0*std::numeric_limits<double>::max();
    for(const auto & val : m_violationField){
        result = std::fmax(result, val);
    }
#if MIMMO_ENABLE_MPI
    MPI_Allreduce(MPI_IN_PLACE, &result, 1, MPI_DOUBLE, MPI_MAX, m_communicator);
#endif
    return    result;
};
 
dmpvector1D *
ControlDeformExtSurface::getViolationField(){
    return  &m_violationField;
};
 
 
double
ControlDeformExtSurface::getTolerance(){
    return m_tolerance;
}
 
const ControlDeformExtSurface::fileListWithType &
ControlDeformExtSurface::getConstraintFiles() const{
    return    m_geoFileList;
}
 
void
ControlDeformExtSurface::setDefField(dmpvecarr3E *field){
    if(!field)  return;
    m_defField.clear();
    m_defField = *field;
};
 
void
ControlDeformExtSurface::setGeometry( mimmo::MimmoSharedPointer<mimmo::MimmoObject> target){
    if(target == nullptr)       return;
    if(target->getType() != 1){
        (*m_log)<<"warning: ControlDeformExtSurface cannot support current geometry. It works only w/ 3D surface."<<std::endl;
        return;
    }
    m_geometry = target;
 
};
 
void
ControlDeformExtSurface::setTolerance(double tol){
    m_tolerance = std::max(0.0, tol);
}
 
void
ControlDeformExtSurface::addConstraint( mimmo::MimmoSharedPointer<mimmo::MimmoObject> constraint){
    if(constraint == nullptr)       return;
    if(constraint->getType() != 1){
        (*m_log)<<"warning: ControlDeformExtSurface cannot support current constraint geometry. It works only w/ 3D surface."<<std::endl;
        return;
    }
    m_geoList.insert(constraint);
};
 
 
 
void
ControlDeformExtSurface::setConstraintFiles(ControlDeformExtSurface::fileListWithType files){
    for(auto & val : files){
        addConstraintFile(val.first, val.second);
    }
};
 
void
ControlDeformExtSurface::addConstraintFile(std::string file, int format){
    if(m_allowed.count(format)>0){
        m_geoFileList[file] = format;
    }
};
 
void     ControlDeformExtSurface::addConstraintFile(std::string file, FileType format){
    addConstraintFile(file,format._to_integral());
};
 
void
ControlDeformExtSurface::removeConstraintFile(std::string file){
    if(m_geoFileList.count(file) >0)    m_geoFileList.erase(file);
};
 
void
ControlDeformExtSurface::removeConstraintFiles(){
    m_geoFileList.clear();
};
 
void
ControlDeformExtSurface::clear(){
    removeConstraintFiles();
    m_geoList.clear();
    m_defField.clear();
    m_violationField.clear();
    m_tolerance = 0.0;
    BaseManipulation::clear();
};
 
void
ControlDeformExtSurface::execute(){
 
    MimmoSharedPointer<MimmoObject> geo = getGeometry();
    if(geo == nullptr){
        (*m_log)<<"Error "+ m_name + " : null pointer to linked geometry found"<<std::endl;
        throw std::runtime_error("Error "+m_name + " : null pointer to linked geometry found");
 
    }
    if(geo->isEmpty()){
        (*m_log)<<"Warning in " + m_name +" : empty linked geometry found"<<std::endl;
    }
 
    bool check = m_defField.getGeometry() == geo;
    check = check && m_defField.getDataLocation() == MPVLocation::POINT;
    if(!check){
        (*m_log)<<"Warning in "<<m_name<<": Unsuitable deformation field linked"<<std::endl;
    }
    m_defField.completeMissingData({{0.0,0.0,0.0}});
 
    m_violationField.clear();
    m_violationField.initialize(geo, MPVLocation::POINT, -1.0E25);
 
    std::vector<long> pointIds = geo->getVerticesIds();
    dvecarr3E points(pointIds.size());
 
    //adding deformation to points and calculate their local
    //barycenter
    int count = 0;
    darray3E dbMin, dbMax;
    dbMin.fill(std::numeric_limits<double>::max());
    dbMax.fill(-1.0*std::numeric_limits<double>::max());
    for (long id : pointIds){
        pointIds[count] = id;
        points[count] = geo->getVertexCoords(id) + m_defField[id];
        for (int i=0; i<3; ++i){
            dbMin[i] = std::min(dbMin[i], points[count][i]);
            dbMax[i] = std::max(dbMax[i], points[count][i]);
        }
        ++count;
    }
    darray3E loc_dcenter = 0.5*(dbMin + dbMax);
 
 
    //evaluate global AABB of geometry and save its center
    darray3E pbMin, pbMax;
    getGlobalBoundingBox(geo, pbMin, pbMax);
    darray3E originalGlobalCenter = 0.5*(pbMin + pbMax);
 
    //put together constrained surfaces in a unique list
    std::vector<MimmoSharedPointer<MimmoObject>> constraint_geos;
    //from files first
    readFileConstraints(constraint_geos);
    // then add those directly linked with addConstraint method.
    constraint_geos.insert(constraint_geos.end(), m_geoList.begin(), m_geoList.end());
 
    //check list size of constraints and rise a warning if the list is empty.
    if(constraint_geos.empty()) {
        (*m_log)<<"Warning in " + m_name +" : no valid constraint geometries are linked to the class. "<<std::endl;
    }
 
    // start examining one by one all external constraints
    darray3E bbMin, bbMax;
 
    for(MimmoSharedPointer<MimmoObject> localg : constraint_geos){
 
        //check constraint skdtree
        localg->buildSkdTree();
        //get the global bounding box of the constraint surface
        getGlobalBoundingBox(localg, bbMin, bbMax);
        //evaluate a guess of the search radius as distance using AABBs of the global
        //target and the global constraint.
        darray3E insMin, insMax;
        double suppval(1.0E-08);
        if(bitpit::CGElem::intersectBoxBox(pbMin,pbMax, bbMin, bbMax, insMin, insMax)){
           suppval =  std::max(1.0E-08, 0.5*norm2(insMax - insMin) );
        }
        double searchRadius = std::max( suppval, norm2(originalGlobalCenter - 0.5*(bbMin + bbMax)) );
 
        //calculate where the center of the global target is located w.r.t to constraint, to keep up the
        //right sign for the distance.
        double refsign = 1.0;
        std::vector<double> distances;
        distances.resize(1);
        evaluateSignedDistance(std::vector<darray3E>(1,originalGlobalCenter), localg, searchRadius, distances);
        if(distances[0] > 0.0) refsign *= -1.0;
 
        //Calculate now the real distances of the deformed points set.
        distances.resize(points.size());
        //evaluate a first guess of the search radius as distance using AABBs of the local
        //point set and the global constraint.
        suppval = 1.0E-08;
        if(bitpit::CGElem::intersectBoxBox(dbMin,dbMax, bbMin, bbMax, insMin, insMax)){
           suppval =  std::max(1.0E-08, 0.5*norm2(insMax - insMin) );
        }
        searchRadius = std::max( suppval, norm2(loc_dcenter - 0.5*(bbMin + bbMax)) );
 
        evaluateSignedDistance(points, localg, searchRadius, distances);
 
        distances *= refsign;
 
        //push values directly in m_violationField.
        count = 0;
        for(long id : pointIds){
            m_violationField[id] = std::max( m_violationField[id], (distances[count] + m_tolerance) );
            ++count;
        }
 
    }//end looping on constraint geometries.
 
    //write resume file.
    writeLog();
 
};
 
void
ControlDeformExtSurface::absorbSectionXML(const bitpit::Config::Section & slotXML, std::string name){
 
    BITPIT_UNUSED(name);
 
    BaseManipulation::absorbSectionXML(slotXML, name);
 
    std::unordered_map<std::string, int > mapp;
    if(slotXML.hasSection("Files")){
 
        const bitpit::Config::Section & filesXML = slotXML.getSection("Files");
 
        for(auto & subfile : filesXML.getSections()){
            std::string path;
            std::string tag, tolstring;
 
            if(subfile.second->hasOption("fullpath"))    {
                path = subfile.second->get("fullpath");
                path = bitpit::utils::string::trim(path);
            }
            if(subfile.second->hasOption("tag")){
                tag = subfile.second->get("tag");
                tag = bitpit::utils::string::trim(tag);
                //check tag;
                auto maybe_tag = FileType::_from_string_nothrow(tag.c_str());
                if(!maybe_tag)    tag.clear();
                else    tag = maybe_tag->_to_string();
            }
            if(!path.empty() && !tag.empty()){
                mapp[path] = int(FileType::_from_string(tag.c_str()));
            }
        }
 
        setConstraintFiles(mapp);
 
    }
 
    if(slotXML.hasOption("Tolerance")){
        std::string input = slotXML.get("Tolerance");
        input = bitpit::utils::string::trim(input);
        double value = 0.0;
        if(!input.empty()){
            std::stringstream ss(input);
            ss >> value;
        }
        setTolerance(value);
    }
 
};
 
void
ControlDeformExtSurface::flushSectionXML(bitpit::Config::Section & slotXML, std::string name){
 
    BITPIT_UNUSED(name);
 
    BaseManipulation::flushSectionXML(slotXML, name);
    bitpit::Config::Section & filesXML = slotXML.addSection("Files");
 
    int counter = 0;
    for(auto & file : m_geoFileList){
        std::string name = "file"+std::to_string(counter);
        bitpit::Config::Section & local = filesXML.addSection(name);
        local.set("fullpath", file.first);
        std::string typetag = (FileType::_from_integral(file.second))._to_string();
        local.set("tag", typetag);
        ++counter;
    }
 
    slotXML.set("Tolerance", std::to_string(m_tolerance));
 
};
 
void
ControlDeformExtSurface::readFileConstraints(std::vector<MimmoSharedPointer<MimmoObject> > & extFileGeos){
 
    extFileGeos.resize(m_geoFileList.size());
 
    int counter = 0;
    std::unique_ptr<MimmoGeometry> geo (new MimmoGeometry(MimmoGeometry::IOMode::READ));
 
    for(auto & geoinfo : m_geoFileList){
 
        svector1D info = extractInfo(geoinfo.first);
        geo->setDir(info[0]);
        geo->setFilename(info[1]);
        geo->setFileType(geoinfo.second);
        geo->execute();
 
        MimmoSharedPointer<MimmoObject> locMesh = geo->getGeometry();
        if (locMesh == nullptr) continue;
 
        extFileGeos[counter] = locMesh;
 
        extFileGeos[counter]->updateAdjacencies();
        extFileGeos[counter]->buildSkdTree();
 
        ++counter;
    }
    extFileGeos.resize(counter);
};
 
svector1D
ControlDeformExtSurface::extractInfo(std::string file){
 
    std::string root, name, tag,temp;
    std::string key1=".", key2="/\\";
 
    std::size_t found = file.find_last_of(key2);
    root = file.substr(0, found);
    temp = file.substr(found+1);
 
    found = temp.find_last_of(key1);
    name = temp.substr(0,found);
    tag = temp.substr(found+1);
 
    svector1D result(3);
    result[0] = root;
    result[1] = name;
    result[2] = tag;
 
    return     result;
}
 
void
ControlDeformExtSurface::evaluateSignedDistance(const std::vector<darray3E> &points, MimmoSharedPointer<MimmoObject> &geo,
                                                double initRadius, std::vector<double> & distances)
{
 
    geo->buildSkdTree();
 
    double rate = 0.05;
    int kmax = 1000;
    int kiter = 0;
    double sRadius(initRadius);
    distances.resize(points.size());
 
    std::vector<darray3E> work = points;
    std::unordered_map<std::size_t,std::size_t> mapPosIndex;
    for(std::size_t i=0; i<points.size(); ++i){
        mapPosIndex[i] = i;
    }
 
    bool checkEmptyWork = work.empty();
#if MIMMO_ENABLE_MPI
    MPI_Allreduce(MPI_IN_PLACE, &checkEmptyWork, 1, MPI_C_BOOL, MPI_LAND, m_communicator);
#endif
 
    while(!checkEmptyWork && kiter < kmax){
 
        std::vector<long> suppCellIds(work.size());
        std::vector<double> distanceWork(work.size());
        std::vector<darray3E> normals(work.size());
 
#if MIMMO_ENABLE_MPI
        std::vector<int> suppCellRanks(work.size());
        skdTreeUtils::signedGlobalDistance(work.size(), work.data(), geo->getSkdTree(), suppCellIds.data(), suppCellRanks.data(), normals.data(), distanceWork.data(), sRadius, false);
#else
        skdTreeUtils::signedDistance(work.size(), work.data(), geo->getSkdTree(), suppCellIds.data(), distanceWork.data(), normals.data(), sRadius);
#endif
 
        //get all points with distances not calculated.
        std::vector<std::array<double,3>> failedPoints;
        std::unordered_map<std::size_t, std::size_t> failedPosIndex;
        failedPoints.reserve(work.size());
        int count(0);
        for(long & val : suppCellIds){
            if(val == bitpit::Cell::NULL_ID){
                failedPoints.push_back(work[count]);
                failedPosIndex[failedPoints.size()-1] = mapPosIndex[count];
            }else{
                distances[mapPosIndex[count]] = distanceWork[count];
            }
            ++count;
        }
 
        std::swap(failedPoints, work);
        std::swap(failedPosIndex, mapPosIndex);
 
        //check work if empty
        checkEmptyWork = work.empty();
#if MIMMO_ENABLE_MPI
        MPI_Allreduce(MPI_IN_PLACE, &checkEmptyWork, 1, MPI_C_BOOL, MPI_LAND, m_communicator);
#endif
        //increase search radius
        sRadius *= (1.0 + rate);
        //increase iteration
        kiter++;
    }
 
    //check for residual unchecked distances, and put them to -1.0E+25
    // just to be manageable by visualizers like paraview.
    for(auto & tuple : mapPosIndex){
        distances[tuple.second] = -1.0E+25;
    }
 
}
 
void
ControlDeformExtSurface::plotOptionalResults(){
 
    m_defField.setName("Deformation Field");
    m_violationField.setName("Violation Distance Field");
    write(getGeometry(), m_defField, m_violationField);
 
}
 
void
ControlDeformExtSurface::writeLog(){
 
    double violationMax = getViolation();
    std::size_t geoListSize = m_geoList.size();
 
#if MIMMO_ENABLE_MPI
    if(m_rank == 0){
#endif
        std::string logname = m_name+std::to_string(getId())+"_violation.log";
        std::ofstream log;
        log.open(logname);
        log<<"mimmo "<<m_name<<" resume file"<<std::endl;
        log<<std::endl;
        log<<std::endl;
        log<<" violation value : " << violationMax << std::endl;
        log<<std::endl;
        log<<" directly linked constraints (number) : " << geoListSize << std::endl;
        log<<std::endl;
        log<<" additional constraints from files : " ;
        for (auto ig : m_geoFileList){
            log<< ig.first <<"  "<<std::endl;
        }
        log.close();
 
#if MIMMO_ENABLE_MPI
    }
    MPI_Barrier(m_communicator); // other ranks stalled here, waiting 0 to finish writing.
 #endif
 
}
 
 
void
ControlDeformExtSurface::getGlobalBoundingBox(MimmoSharedPointer<MimmoObject> & geo, darray3E & bMin, darray3E & bMax){
    if(geo->getBoundingBoxSyncStatus() != mimmo::SyncStatus::SYNC){
        geo->update();
    }
    geo->getBoundingBox(bMin, bMax, true);
}
 
}//end mimmo namespace