TwistGeometry.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 "TwistGeometry.hpp"
 
namespace mimmo{
 
 
TwistGeometry::TwistGeometry(darray3E origin, darray3E direction){
    m_origin = origin;
    m_direction = direction;
    m_sym = false;
    m_alpha = 0.0;
    m_distance = 0.0;
    m_name = "mimmo.TwistGeometry";
};
 
TwistGeometry::TwistGeometry(const bitpit::Config::Section & rootXML){
 
    m_origin.fill(0.0);
    m_direction.fill(0.0);
    m_alpha = 0.0;
    m_distance = 0.0;
    m_sym = false;
    m_name = "mimmo.TwistGeometry";
 
    std::string fallback_name = "ClassNONE";
    std::string input = rootXML.get("ClassName", fallback_name);
    input = bitpit::utils::string::trim(input);
    if(input == "mimmo.TwistGeometry"){
        absorbSectionXML(rootXML);
    }else{
        warningXML(m_log, m_name);
    };
}
 
TwistGeometry::~TwistGeometry(){};
 
TwistGeometry::TwistGeometry(const TwistGeometry & other):BaseManipulation(other){
    m_origin = other.m_origin;
    m_direction = other.m_direction;
    m_alpha = other.m_alpha;
    m_distance = other.m_distance;
    m_sym = other.m_sym;
    m_filter = other.m_filter;
};
 
TwistGeometry& TwistGeometry::operator=(TwistGeometry other){
    swap(other);
    return *this;
};
 
void TwistGeometry::swap(TwistGeometry & x) noexcept
{
    std::swap(m_origin, x.m_origin);
    std::swap(m_direction, x.m_direction);
    std::swap(m_alpha, x.m_alpha);
    std::swap(m_distance, x.m_distance);
    std::swap(m_sym, x.m_sym);
    m_filter.swap(x.m_filter);
    m_displ.swap(x.m_displ);
    BaseManipulation::swap(x);
};
 
void
TwistGeometry::buildPorts(){
    bool built = true;
    built = (built && createPortIn<darray3E, TwistGeometry>(&m_origin, M_POINT));
    built = (built && createPortIn<darray3E, TwistGeometry>(&m_direction, M_AXIS));
    built = (built && createPortIn<double, TwistGeometry>(&m_alpha, M_VALUED));
    built = (built && createPortIn<double, TwistGeometry>(&m_distance, M_VALUED2));
    built = (built && createPortIn<dmpvector1D*, TwistGeometry>(this, &mimmo::TwistGeometry::setFilter, M_FILTER));
    built = (built && createPortIn<MimmoSharedPointer<MimmoObject>, TwistGeometry>(&m_geometry, M_GEOM, true));
    built = (built && createPortOut<dmpvecarr3E*, TwistGeometry>(this, &mimmo::TwistGeometry::getDisplacements, M_GDISPLS));
    built = (built && createPortOut<MimmoSharedPointer<MimmoObject>, TwistGeometry>(this, &BaseManipulation::getGeometry, M_GEOM));
    m_arePortsBuilt = built;
};
 
void
TwistGeometry::setAxis(darray3E origin, darray3E direction){
    m_origin = origin;
    m_direction = direction;
}
 
void
TwistGeometry::setOrigin(darray3E origin){
    m_origin = origin;
}
 
void
TwistGeometry::setDirection(darray3E direction){
    m_direction = direction;
    double L = norm2(m_direction);
    for (int i=0; i<3; i++)
        m_direction[i] /= L;
}
 
void
TwistGeometry::setTwist(double alpha){
    m_alpha = alpha;
}
 
void
TwistGeometry::setSym(bool sym){
    m_sym = sym;
}
 
void
TwistGeometry::setMaxDistance(double distance){
    m_distance = distance;
}
 
void
TwistGeometry::setFilter(dmpvector1D *filter){
    if(!filter) return;
    m_filter = *filter;
}
 
dmpvecarr3E*
TwistGeometry::getDisplacements(){
    return &m_displ;
};
 
void
TwistGeometry::execute(){
 
    if(getGeometry() == nullptr){
        (*m_log)<<m_name + " : nullptr pointer to linked geometry found"<<std::endl;
        throw std::runtime_error(m_name + "nullptr pointer to linked geometry found");
    }
 
    checkFilter();
 
    m_displ.clear();
    m_displ.setDataLocation(mimmo::MPVLocation::POINT);
    m_displ.reserve(getGeometry()->getNVertices());
    m_displ.setGeometry(getGeometry());
 
    darray3E point, rotated;
    long ID;
    darray3E projected;
    double distance;
    //double sign;
    double rot;
    darray3E value;
 
    for (const auto & vertex : m_geometry->getVertices()){
        point = vertex.getCoords();
        ID = vertex.getId();
 
        //signed distance from origin
        distance = dotProduct((point-m_origin),m_direction);
 
        //compute coefficients and constant vectors of rodriguez formula
        rot = std::min(m_alpha, (std::abs(distance)/m_distance)*m_alpha);
        if (distance < 0) rot = -rot*int(m_sym);
        double a = cos(rot);
        darray3E b =  (1 - cos(rot)) * m_direction;
        double c = sin(rot);
 
        //project point on axis (local origin)
        projected = distance*m_direction + m_origin;
 
 
        point -= projected;
        //rodrigues formula
        rotated = a * point +
                b * dotProduct(m_direction, point) +
                c * crossProduct(m_direction, point);
 
        rotated += projected;
        point += projected;
 
        value = (rotated-point)*m_filter[ID];
        m_displ.insert(ID, value);
 
    }
};
 
void
TwistGeometry::apply(){
    _apply(m_displ);
}
 
void
TwistGeometry::checkFilter(){
 
    bool check = m_filter.getDataLocation() == mimmo::MPVLocation::POINT;
    check = check && m_filter.completeMissingData(0.0);
    check = check && m_filter.getGeometry() == getGeometry();
 
    if (!check){
        m_log->setPriority(bitpit::log::Verbosity::DEBUG);
        (*m_log)<<"Not valid filter found in "<<m_name<<". Proceeding with default unitary field"<<std::endl;
        m_log->setPriority(bitpit::log::Verbosity::NORMAL);
 
        m_filter.clear();
        m_filter.setGeometry(m_geometry);
        m_filter.setDataLocation(mimmo::MPVLocation::POINT);
        m_filter.reserve(getGeometry()->getNVertices());
        for (const auto & vertex : getGeometry()->getVertices()){
            m_filter.insert(vertex.getId(), 1.0);
        }
    }
}
 
void
TwistGeometry::absorbSectionXML(const bitpit::Config::Section & slotXML, std::string name){
 
    BITPIT_UNUSED(name);
 
    BaseManipulation::absorbSectionXML(slotXML, name);
 
    if(slotXML.hasOption("Origin")){
        std::string input = slotXML.get("Origin");
        input = bitpit::utils::string::trim(input);
        darray3E temp = {{0.0,0.0,0.0}};
        if(!input.empty()){
            std::stringstream ss(input);
            for(auto &val : temp) ss>>val;
        }
        setOrigin(temp);
    }
 
    if(slotXML.hasOption("Direction")){
        std::string input = slotXML.get("Direction");
        input = bitpit::utils::string::trim(input);
        darray3E temp = {{0.0,0.0,0.0}};
        if(!input.empty()){
            std::stringstream ss(input);
            for(auto &val : temp) ss>>val;
        }
        setDirection(temp);
    }
 
    if(slotXML.hasOption("Twist")){
        std::string input = slotXML.get("Twist");
        input = bitpit::utils::string::trim(input);
        double temp = 0.0;
        if(!input.empty()){
            std::stringstream ss(input);
            ss>>temp;
        }
        setTwist(temp);
    }
 
    if(slotXML.hasOption("Distance")){
        std::string input = slotXML.get("Distance");
        input = bitpit::utils::string::trim(input);
        double temp = 0.0;
        if(!input.empty()){
            std::stringstream ss(input);
            ss>>temp;
        }
        setMaxDistance(temp);
    }
 
    if(slotXML.hasOption("Symmetric")){
        std::string input = slotXML.get("Symmetric");
        input = bitpit::utils::string::trim(input);
        bool temp = false;
        if(!input.empty()){
            std::stringstream ss(input);
            ss>>temp;
        }
        setSym(temp);
    };
 
};
 
void
TwistGeometry::flushSectionXML(bitpit::Config::Section & slotXML, std::string name){
 
    BITPIT_UNUSED(name);
 
    BaseManipulation::flushSectionXML(slotXML, name);
 
    {
        std::stringstream ss;
        ss<<std::scientific<<m_origin[0]<<'\t'<<m_origin[1]<<'\t'<<m_origin[2];
        slotXML.set("Origin", ss.str());
    }
 
    {
        std::stringstream ss;
        ss<<std::scientific<<m_direction[0]<<'\t'<<m_direction[1]<<'\t'<<m_direction[2];
        slotXML.set("Direction", ss.str());
    }
 
    slotXML.set("Twist", std::to_string(m_alpha));
 
    slotXML.set("Distance", std::to_string(m_distance));
 
    slotXML.set("Symmetric", std::to_string(int(m_sym)));
 
};
 
}