ClipGeometry.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 "ClipGeometry.hpp"
 
namespace mimmo{
 
ClipGeometry::ClipGeometry(){
    m_name = "mimmo.ClipGeometry";
    m_plane.fill(0.0);
    m_origin.fill(0.0);
    m_normal.fill(0.0);
    m_insideout = false;
    m_clipped.reset();
    m_implicit = false;
};
 
ClipGeometry::ClipGeometry(const bitpit::Config::Section & rootXML){
 
    m_name = "mimmo.ClipGeometry";
    m_plane.fill(0.0);
    m_insideout = false;
    m_clipped.reset();
    m_origin.fill(0.0);
    m_normal.fill(0.0);
    m_implicit = false;
 
    std::string fallback_name = "ClassNONE";
    std::string input = rootXML.get("ClassName", fallback_name);
    input = bitpit::utils::string::trim(input);
    if(input == "mimmo.ClipGeometry"){
        absorbSectionXML(rootXML);
    }else{
        warningXML(m_log, m_name);
    };
}
 
ClipGeometry::~ClipGeometry(){};
 
ClipGeometry::ClipGeometry(const ClipGeometry & other):BaseManipulation(other){
    m_plane = other.m_plane;
    m_insideout = other.m_insideout;
    m_origin = other.m_origin;
    m_normal = other.m_normal;
    m_implicit = other.m_implicit;
};
 
void
ClipGeometry::buildPorts(){
 
    bool built = true;
 
    built = (built && createPortIn<darray4E, ClipGeometry>(this, &mimmo::ClipGeometry::setClipPlane, M_PLANE));
    built = (built && createPortIn<darray3E, ClipGeometry>(this, &mimmo::ClipGeometry::setOrigin, M_POINT));
    built = (built && createPortIn<darray3E, ClipGeometry>(this, &mimmo::ClipGeometry::setNormal, M_AXIS));
    built = (built && createPortIn<bool, ClipGeometry>(this, &mimmo::ClipGeometry::setInsideOut, M_VALUEB));
    built = (built && createPortIn<MimmoSharedPointer<MimmoObject>, ClipGeometry>(this, &mimmo::ClipGeometry::setGeometry, M_GEOM, true));
 
    built = (built && createPortOut<MimmoSharedPointer<MimmoObject>, ClipGeometry>(this, &mimmo::ClipGeometry::getClippedPatch, M_GEOM));
    m_arePortsBuilt = built;
};
 
bool
ClipGeometry::isInsideOut(){
    return(m_insideout);
};
 
MimmoSharedPointer<MimmoObject>
ClipGeometry::getClippedPatch(){
    return(m_clipped);
};
 
darray4E
ClipGeometry::getClipPlane(){
    return m_plane;
};
 
void
ClipGeometry::setClipPlane(darray4E plane){
    m_plane = plane;
    m_implicit = true;
};
 
void
ClipGeometry::setClipPlane(darray3E origin, darray3E normal){
 
    normal /= norm2(normal);
    double b = -1.0*dotProduct(origin, normal);
 
    m_plane[0] = normal[0];
    m_plane[1] = normal[1];
    m_plane[2] = normal[2];
    m_plane[3] = b;
 
    m_implicit = true;
 
};
 
void
ClipGeometry::setOrigin(darray3E origin){
    m_origin = origin;
};
 
void
ClipGeometry::setNormal(darray3E normal){
    m_normal = normal;
};
 
void
ClipGeometry::setInsideOut(bool flag){
    m_insideout = flag;
};
 
void
ClipGeometry::execute(){
 
    if(getGeometry() == nullptr){
        (*m_log)<<m_name + " : nullptr geometry linked."<<std::endl;
        throw std::runtime_error (m_name + " : nullptr geometry linked.");
    };
 
    /* If an implicit definition is not present it has to be computed
     * by using origin and normal.
     */
    if (!m_implicit) setClipPlane(m_origin, m_normal);
 
    m_clipped.reset();
 
    livector1D extracted = clipPlane();
 
    /* Create subpatch.*/
    MimmoSharedPointer<MimmoObject> temp(new MimmoObject(getGeometry()->getType()));
    m_clipped = temp;
    bitpit::PatchKernel * tri = getGeometry()->getPatch();
 
    livector1D idVertexList = getGeometry()->getVertexFromCellList(extracted);
 
    temp->getPatch()->reserveVertices(idVertexList.size());
    temp->getPatch()->reserveCells(extracted.size());
 
    for(long idV : idVertexList){
        m_clipped->addVertex(tri->getVertexCoords(idV),idV);
    }
 
    int rank;
    for(long idCell : extracted){
 
        bitpit::Cell & cell = tri->getCell(idCell);
        rank  =-1;
#if MIMMO_ENABLE_MPI
        rank = getGeometry()->getPatch()->getCellRank(idCell);
#endif
        m_clipped->addCell(cell, idCell, rank);
    }
    //MPI version: recovering info from the local rank partition
    auto originalmap = getGeometry()->getPIDTypeListWNames();
    auto currentPIDmap = temp->getPIDTypeList();
    for(const auto & val: currentPIDmap){
        m_clipped->setPIDName(val, originalmap[val]);
    }
 
    // Clean and Update mesh. This will update even parallel structures if needed
    m_clipped->cleanGeometry();
    m_clipped->update();
 
};
 
livector1D
ClipGeometry::clipPlane(){
 
    livector1D result;
    darray3E norm;
    double offset;
    int counter;
    double sig = 1.0 - 2.0*(int)isInsideOut();
    long iD;
 
 
    for(int i=0; i<3; ++i)norm[i] = m_plane[i];
    offset = m_plane[3];
 
    double normPlane = norm2(norm);
    if(normPlane < 1.E-18)return result;
    norm /= normPlane;
 
    bitpit::PatchKernel * tri = getGeometry()->getPatch();
 
    if(getGeometry()->getType() == 3){
        counter = 0;
        result.resize(tri->getVertexCount());
        for(const auto vert : tri->getVertices()){
            iD = vert.getId();
            if(sig*(dotProduct(norm, vert.getCoords()) + offset) >0){
                result[counter] = iD;
                ++counter;
            }
        }
 
    }else{
        counter = 0;
        result.resize(tri->getCellCount());
        for(const auto cell : tri->getCells()){
            iD = cell.getId();
            if(sig*(dotProduct(norm, tri->evalCellCentroid(iD)) + offset) >0){
                result[counter] = iD;
                ++counter;
            }
        }
    }
    result.resize(counter);
    return result;
};
 
void
ClipGeometry::plotOptionalResults(){
 
    write(getClippedPatch());
 
}
 
void
ClipGeometry::absorbSectionXML(const bitpit::Config::Section & slotXML, std::string name){
 
    /*start absorbing*/
    BITPIT_UNUSED(name);
 
    BaseManipulation::absorbSectionXML(slotXML,name);
 
    if(slotXML.hasOption("InsideOut")){
        std::string input = slotXML.get("InsideOut");
        input = bitpit::utils::string::trim(input);
        bool value = false;
        if(!input.empty()){
            std::stringstream ss(input);
            ss >> value;
        }
        setInsideOut(value);
    }
 
    if(slotXML.hasSection("ClipPlane")){
        const bitpit::Config::Section & planeXML = slotXML.getSection("ClipPlane");
 
        std::string input1 = planeXML.get("Point");
        std::string input2 = planeXML.get("Normal");
        input1 = bitpit::utils::string::trim(input1);
        input2 = bitpit::utils::string::trim(input2);
 
        darray3E temp1 = {{0.0,0.0,0.0}};
        darray3E temp2 = {{0.0,0.0,0.0}};
 
        if(!input1.empty()){
            std::stringstream ss(input1);
            ss>>temp1[0]>>temp1[1]>>temp1[2];
        }
        if(!input2.empty()){
            std::stringstream ss(input2);
            ss>>temp2[0]>>temp2[1]>>temp2[2];
        }
 
        setClipPlane(temp1, temp2);
    };
};
 
void
ClipGeometry::flushSectionXML(bitpit::Config::Section & slotXML, std::string name){
 
    BITPIT_UNUSED(name);
 
    BaseManipulation::flushSectionXML(slotXML, name);
 
    int value = m_insideout;
    slotXML.set("InsideOut", std::to_string(value));
 
    {
        darray4E org = getClipPlane();
        darray3E normal;
        darray3E point = {{0.0,0.0,0.0}};
        int imax = -1;
        double dum = 0.0;
        for(int i=0; i<3; ++i){
            normal[i] =org[i];
            if(abs(normal[i]) > dum) {
                imax = i;
                dum = abs(normal[i]);
            }
        }
        if(imax != -1)point[imax] = -1.0*org[3]/normal[imax];
 
        std::stringstream ss1, ss2;
        ss1<<std::scientific<<point[0]<<'\t'<<point[1]<<'\t'<<point[2];
        ss2<<std::scientific<<normal[0]<<'\t'<<normal[1]<<'\t'<<normal[2];
 
        bitpit::Config::Section & planeXML = slotXML.addSection("ClipPlane");
        planeXML.set("Point",ss1.str());
        planeXML.set("Normal",ss2.str());
    }
 
};
 
}