SpecularPoints.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 "SpecularPoints.hpp"
#include "SkdTreeUtils.hpp"
 
namespace mimmo{
 
SpecularPoints::SpecularPoints(){
    m_name = "mimmo.SpecularPoints";
    m_plane.fill(0.0);
    m_insideout = false;
    m_force = true;
    m_origin.fill(0.0);
    m_normal.fill(0.0);
    m_implicit = false;
    m_scalar = nullptr;
    m_vector = nullptr;
};
 
SpecularPoints::SpecularPoints(const bitpit::Config::Section & rootXML){
 
    m_name = "mimmo.SpecularPoints";
    m_plane.fill(0.0);
    m_insideout = false;
    m_force = true;
    m_origin.fill(0.0);
    m_normal.fill(0.0);
    m_implicit = false;
    m_scalar = nullptr;
    m_vector = nullptr;
 
    std::string fallback_name = "ClassNONE";
    std::string input = rootXML.get("ClassName", fallback_name);
    input = bitpit::utils::string::trim(input);
    if(input == "mimmo.SpecularPoints"){
        absorbSectionXML(rootXML);
    }else{
        warningXML(m_log, m_name);
    };
}
 
SpecularPoints::~SpecularPoints(){};
 
SpecularPoints::SpecularPoints(const SpecularPoints & other):ProjPatchOnSurface(other){
    m_plane = other.m_plane;
    m_scalar = other.m_scalar;
    m_vector = other.m_vector;
    m_insideout = other.m_insideout;
    m_force = other.m_force;
    m_origin = other.m_origin;
    m_normal = other.m_normal;
    m_implicit = other.m_implicit;
    m_pc = other.m_pc;
};
 
SpecularPoints & SpecularPoints::operator=(SpecularPoints other){
    swap(other);
    return *this;
}
 
void SpecularPoints::swap(SpecularPoints & x) noexcept
{
    std::swap(m_plane, x.m_plane);
    std::swap(m_scalar, x.m_scalar);
    std::swap(m_vector, x.m_vector);
    std::swap(m_insideout, x.m_insideout);
    std::swap(m_force, x.m_force);
    std::swap(m_origin, x.m_origin);
    std::swap(m_normal, x.m_normal);
    std::swap(m_implicit, x.m_implicit);
    m_scalarMirrored.swap(x.m_scalarMirrored);
    m_vectorMirrored.swap(x.m_vectorMirrored);
    std::swap(m_pc,x.m_pc);
 
    ProjPatchOnSurface::swap(x);
}
void
SpecularPoints::buildPorts(){
    bool built = true;
 
    built = (built && createPortIn<MimmoSharedPointer<MimmoObject>, SpecularPoints>(this,&mimmo::SpecularPoints::setPointCloud, M_GEOM2, true));
    built = (built && createPortIn<dmpvecarr3E*, SpecularPoints>(this, &mimmo::SpecularPoints::setVectorData, M_VECTORFIELD));
    built = (built && createPortIn<dmpvector1D*, SpecularPoints>(this, &mimmo::SpecularPoints::setScalarData, M_DATAFIELD));
    built = (built && createPortIn<darray4E, SpecularPoints>(this, &mimmo::SpecularPoints::setPlane, M_PLANE));
    built = (built && createPortIn<darray3E, SpecularPoints>(this, &mimmo::SpecularPoints::setOrigin, M_POINT));
    built = (built && createPortIn<darray3E, SpecularPoints>(this, &mimmo::SpecularPoints::setNormal, M_AXIS));
    built = (built && createPortIn<MimmoSharedPointer<MimmoObject>, SpecularPoints>(this, &mimmo::SpecularPoints::setGeometry, M_GEOM));
 
    built = (built && createPortOut<MimmoSharedPointer<MimmoObject>, SpecularPoints>(this, &mimmo::SpecularPoints::getMirroredPointCloud, M_GEOM));
    built = (built && createPortOut<dmpvecarr3E*, SpecularPoints>(this, &mimmo::SpecularPoints::getMirroredVectorData, M_VECTORFIELD));
    built = (built && createPortOut<dmpvector1D*, SpecularPoints>(this, &mimmo::SpecularPoints::getMirroredScalarData, M_SCALARFIELD));
    m_arePortsBuilt = built;
};
 
dmpvector1D *
SpecularPoints::getOriginalScalarData(){
    return m_scalar;
};
 
dmpvecarr3E *
SpecularPoints::getOriginalVectorData(){
    return m_vector;
};
 
dmpvector1D *
SpecularPoints::getMirroredScalarData(){
    return &m_scalarMirrored;
};
 
dmpvecarr3E *
SpecularPoints::getMirroredVectorData(){
    return &m_vectorMirrored;
};
 
 
dvecarr3E
SpecularPoints::getMirroredRawCoords(){
 
    dvecarr3E points;
    if(m_patch != nullptr){
        points.reserve(m_patch->getNVertices());
        for(bitpit::Vertex & vert : m_patch->getVertices()){
            points.push_back(vert.getCoords());
        }
    }
    return points;
};
 
livector1D
SpecularPoints::getMirroredLabels(){
 
    livector1D labels;
    if(m_patch != nullptr){
        labels.reserve(m_patch->getNVertices());
        for(bitpit::Vertex & vert : m_patch->getVertices()){
            labels.push_back(vert.getId());
        }
    }
    return labels;
};
 
MimmoSharedPointer<MimmoObject>
SpecularPoints::getMirroredPointCloud(){
    return m_patch;
}
 
darray4E
SpecularPoints::getPlane(){
    return  m_plane;
};
 
bool
SpecularPoints::isInsideOut(){
    return m_insideout;
}
 
bool
SpecularPoints::isForced(){
    return m_force;
}
 
 
void
SpecularPoints::setPointCloud(MimmoSharedPointer<MimmoObject> targetpatch){
 
    if(targetpatch == nullptr) return;
    if(targetpatch->getType() !=3) return;
    m_pc = targetpatch;
};
 
 
void
SpecularPoints::setScalarData(dmpvector1D * scalardata){
    if(scalardata == nullptr) return;
    if(scalardata->getDataLocation() != MPVLocation::POINT)  return;
    m_scalar = scalardata;
};
 
void
SpecularPoints::setVectorData(dmpvecarr3E * vectordata){
    if(vectordata == nullptr) return;
    if(vectordata->getDataLocation() != MPVLocation::POINT)  return;
    m_vector = vectordata;
};
 
void
SpecularPoints::setPlane(darray4E plane){
    m_plane = plane;
    m_implicit = true;
};
 
void
SpecularPoints::setPlane(darray3E origin, darray3E normal){
 
    double normx=norm2(normal);
    if(normx > std::numeric_limits<double>::min())  normal /= normx;
 
    m_plane[0] = normal[0];
    m_plane[1] = normal[1];
    m_plane[2] = normal[2];
    m_plane[3] = -1.0*dotProduct(origin, normal);
 
    m_implicit = true;
 
};
 
void
SpecularPoints::setOrigin(darray3E origin){
    m_origin = origin;
};
 
void
SpecularPoints::setNormal(darray3E normal){
    m_normal = normal;
};
 
void
SpecularPoints::setInsideOut(bool flag){
    m_insideout = flag;
};
 
void
SpecularPoints::setForce(bool flag){
    m_force = flag;
}
 
 
void
SpecularPoints::execute(){
 
    //check if there is a valid a point cloud
    if (m_pc == nullptr){
        (*m_log)<<"Error in "<<m_name << " : nullptr pointer to target point cloud"<<std::endl;
        throw std::runtime_error (m_name + " : nullptr pointer to target point cloud");
    }
 
    //clone it into the internal member m_cobj (mother class ProjPatchOnSurface)
    m_cobj = m_pc->clone();
    //initialize mirrored data;
    m_scalarMirrored.clear();
    m_vectorMirrored.clear();
    m_scalarMirrored.initialize(m_cobj, MPVLocation::POINT, 0.);
    m_scalarMirrored.setName("MirroredScalarData");
    m_vectorMirrored.initialize(m_cobj, MPVLocation::POINT, {{0.0,0.0,0.0}});
    m_vectorMirrored.setName("MirroredVectorData");
 
 
    /* Check the plane stuff and if an implicit definition is not present it has to be computed
     * by using origin and normal.
     */
    if (!m_implicit) setPlane(m_origin, m_normal);
 
    darray3E norm;
    double offsetPlane;
    for(int i=0; i<3; ++i)    norm[i] = m_plane[i];
    offsetPlane = m_plane[3];
    double normPlane = norm2(norm);
 
    //if normal is 0 there is nothing to do
    if(normPlane <= std::numeric_limits<double>::min()){
        (*m_log)<< "Error: "<<m_name <<" : no valid plane normal found"<<std::endl;
        m_patch = m_cobj;
        return;
    }
 
    //adjust the normal to be unitary in modulus
    norm /= normPlane;
 
    //see if need to project the cloud on a surface geometry: set bool project as active.
    bool project = false;
    int cellcount(0);
    if (getGeometry() != nullptr){
        cellcount = getGeometry()->getNCells();
#if MIMMO_ENABLE_MPI
        MPI_Allreduce(MPI_IN_PLACE, &cellcount, 1, MPI_INT, MPI_SUM, m_communicator);
#endif
        project = (cellcount > 0);
    }
 
 
    //calculate a reasonable margin for defining the plane offset. This is
    // done to distinguish quickly points which lies on plane or not.
    double margin = 1.0E-12;
    if(project){
        double areaTot = 0.0;
        bitpit::SurfaceKernel * tri = static_cast<bitpit::SurfaceKernel * >(getGeometry()->getPatch());
        for(auto it = tri->internalCellBegin(); it != tri->internalCellEnd(); ++it ){
            areaTot += tri->evalCellArea(it.getId());
        }
#if MIMMO_ENABLE_MPI
        MPI_Allreduce(MPI_IN_PLACE, &areaTot, 1, MPI_DOUBLE, MPI_SUM, m_communicator);
#endif
        if(areaTot > std::numeric_limits<double>::min()) {
            margin =  1.2*std::pow(areaTot/((double)cellcount),0.5);
        }
    }
 
    //take into account the right emiplane with insideout given by the User.
    double sig = (1.0  - 2.0*((int)m_insideout));
 
    //before mirroring, evaluate starting offsets for labeling new vertices to be added.
 
    livector1D localVertIds = m_cobj->getVerticesIds();
    int nPCvert = int(localVertIds.size());
 
    long offsetVertLabel(0);
    if(!localVertIds.empty())
        offsetVertLabel = *(std::max_element(localVertIds.begin(), localVertIds.end())) + 1;
    long offsetCellLabel(0);
    {
        livector1D localCellIds = m_cobj->getCellsIds();
        if(!localCellIds.empty())
            offsetCellLabel = *(std::max_element(localCellIds.begin(), localCellIds.end())) + 1;
    }
 
 
#if MIMMO_ENABLE_MPI
    std::vector<int> rankPCSize (m_nprocs, 0);
    rankPCSize[m_rank] = nPCvert;
    MPI_Allreduce(MPI_IN_PLACE, &offsetVertLabel, 1, MPI_LONG, MPI_MAX, m_communicator);
    MPI_Allreduce(MPI_IN_PLACE, &offsetCellLabel, 1, MPI_LONG, MPI_MAX, m_communicator);
    MPI_Allreduce(MPI_IN_PLACE, rankPCSize.data(), m_nprocs, MPI_INT, MPI_MAX, m_communicator);
 
    for(int i=0; i<m_rank; ++i) {
        offsetVertLabel += long(rankPCSize[i]);
        offsetCellLabel += long(rankPCSize[i]);
    }
#endif
 
    //copy data of scalar and vector (if any) into mirrored twin structures
    //(that are already initialized to default 0 value) and reserve them
    //for new vertices insertion
    if(m_scalar){
        for(long id: localVertIds){
            if(m_scalar->exists(id)){
                m_scalarMirrored[id] = m_scalar->at(id);
            }
        }
        m_scalarMirrored.reserve(2*nPCvert);
    }
 
    if(m_vector){
        for(long id: localVertIds){
            if(m_vector->exists(id)){
                m_vectorMirrored[id] = m_vector->at(id);
            }
        }
        m_vectorMirrored.reserve(2*nPCvert);
    }
 
    m_cobj->getPatch()->reserveVertices(2*nPCvert);
    m_cobj->getPatch()->reserveCells(2*nPCvert);
 
 
    //FINALLY mirror the points, scalar and vector fields.
    darray3E vert;
    bitpit::ElementType type = bitpit::ElementType::VERTEX;
    for (long id: localVertIds){
        vert = m_cobj->getVertexCoords(id);
        double distance = sig*(dotProduct(norm, vert) + offsetPlane);
 
        if(distance > margin || m_force){
 
            m_cobj->addVertex( (vert - 2.0*distance*sig*norm), offsetVertLabel);
#if MIMMO_ENABLE_MPI
            m_cobj->addConnectedCell(std::vector<long>(1, offsetVertLabel), type, 0, offsetCellLabel, m_rank);
#else
            m_cobj->addConnectedCell(std::vector<long>(1, offsetVertLabel), type, 0, offsetCellLabel);
#endif
            m_scalarMirrored.insert(offsetVertLabel, m_scalarMirrored[id] );
            m_vectorMirrored.insert(offsetVertLabel, (m_vectorMirrored[id] -2.0*dotProduct(m_vectorMirrored[id], sig*norm)*sig*norm) );
 
            ++offsetVertLabel;
            ++offsetCellLabel;
        }
    }
 
    if(project){
        //this chunk will project the original cloud onto the target surface
        //m_cobj is filled up before.
        m_workingOnTarget = true;
        //final m_patch is m_cobj itself;
        projection();
        //for MPI version, this method already update the patch.
    }else{
        m_patch = m_cobj;
        m_patch->cleanPatchInfo();
    #if MIMMO_ENABLE_MPI
        m_patch->resetPointGhostExchangeInfo();
    #endif
        m_patch->update();
    }
 
    //squeeze stuff out
    m_patch->getPatch()->squeezeVertices();
    m_patch->getPatch()->squeezeCells();
    m_scalarMirrored.shrinkToFit();
    m_vectorMirrored.shrinkToFit();
 
    //evaluate kdtree if required by the class
    if(m_patch){
        if(m_buildKdTree)    m_patch->buildKdTree();
    }else{
        (*m_log)<<m_name << " : failed mirroring Point Cloud"<<std::endl;
    }
};
 
void
SpecularPoints::clear(){
    ProjPatchOnSurface::clear();
    m_scalar = nullptr;
    m_vector = nullptr;
    m_scalarMirrored.clear();
    m_vectorMirrored.clear();
    m_plane.fill(0.0);
    m_insideout = false;
};
 
void
SpecularPoints::absorbSectionXML(const bitpit::Config::Section & slotXML, std::string name){
 
    BITPIT_UNUSED(name);
 
    //start absorbing
    ProjPatchOnSurface::absorbSectionXML(slotXML, name);
 
    if(slotXML.hasOption("Force")){
        std::string input = slotXML.get("Force");
        input = bitpit::utils::string::trim(input);
        bool value = false;
        if(!input.empty()){
            std::stringstream ss(input);
            ss >> value;
        }
        setForce(value);
    }
 
    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("Plane")){
        const bitpit::Config::Section & planeXML = slotXML.getSection("Plane");
 
        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];
        }
 
        setPlane(temp1, temp2);
    };
};
 
void
SpecularPoints::flushSectionXML(bitpit::Config::Section & slotXML, std::string name){
 
    BITPIT_UNUSED(name);
 
    ProjPatchOnSurface::flushSectionXML(slotXML, name);
 
    int value = m_force;
    slotXML.set("Force", std::to_string(value));
 
    value = m_insideout;
    slotXML.set("InsideOut", std::to_string(value));
 
 
    {
        darray4E org = getPlane();
        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("Plane");
        planeXML.set("Point",ss1.str());
        planeXML.set("Normal",ss2.str());
    }
 
};
 
void
SpecularPoints::plotOptionalResults(){
    write(m_patch, m_scalarMirrored, m_vectorMirrored);
};
 
}