BasicMeshes.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 "BasicMeshes.hpp"
#include "customOperators.hpp"
#include <Operators.hpp>
#include <bitpit_common.hpp>
 
namespace mimmo{
 
UStructMesh::UStructMesh(){
    m_nx = 0; m_ny=0; m_nz=0;
    m_dx = 0; m_dy=0; m_dz=0;
 
    m_setorigin = false;
    m_setspan = false;
    m_setInfLimits = false;
    m_setRefSys = false;
    m_isBuild = false;
    m_origin_temp = {{0.0,0.0,0.0}};
    m_span_temp = {{1.0,1.0,1.0}};
    m_inflimits_temp = {{0.0,0.0,0.0}};
    m_refsystem_temp[0] = {{1.0,0.0,0.0}};
    m_refsystem_temp[1] = {{0.0,1.0,0.0}};
    m_refsystem_temp[2] = {{0.0,0.0,1.0}};
    m_shapetype_temp = ShapeType::CUBE;
};
 
UStructMesh::~UStructMesh(){}
 
UStructMesh::UStructMesh(const UStructMesh & other){
 
    // Number of cells
    m_nx = other.m_nx;
    m_ny = other.m_ny;
    m_nz = other.m_nz;
 
    // Cell Spacing
    m_dx = other.m_dx;
    m_dy = other.m_dy;
    m_dz = other.m_dz;
 
    // Copy cell edges and cell centers ----------------------------------------- //
    m_xnode = other.m_xnode;
    m_ynode = other.m_ynode;
    m_znode = other.m_znode;
    m_xedge = other.m_xedge;
    m_yedge = other.m_yedge;
    m_zedge = other.m_zedge;
 
    if(other.m_shape){
        switch(other.m_shape->getShapeType()){
            case ShapeType::WEDGE :
            {
                const Wedge * pp = dynamic_cast<const Wedge*>(other.m_shape.get());
                if (pp != nullptr)  m_shape = std::unique_ptr<BasicShape>(new Wedge(*(pp)));
            }
            break;
            case ShapeType::CYLINDER :
                {
                    const Cylinder * pp = dynamic_cast<const Cylinder*>(other.m_shape.get());
                    if (pp != nullptr)  m_shape = std::unique_ptr<BasicShape>(new Cylinder(*(pp)));
                }
                break;
            case ShapeType::SPHERE :
                {
                    const Sphere * pp = dynamic_cast<const Sphere*>(other.m_shape.get());
                    if (pp != nullptr)  m_shape = std::unique_ptr<BasicShape>(new Sphere(*(pp)));
                }
                break;
            default://CUBE
                {
                    const Cube * pp = dynamic_cast<const Cube*>(other.m_shape.get());
                    if (pp != nullptr)  m_shape = std::unique_ptr<BasicShape>(new Cube(*(pp)));
                }
                break;
        }
    }
 
    //copy temporary members
    m_origin_temp = other.m_origin_temp;
    m_span_temp = other.m_span_temp;
    m_inflimits_temp = other.m_inflimits_temp;
    m_refsystem_temp = other.m_refsystem_temp;
    m_shapetype_temp = other.m_shapetype_temp;
    m_setorigin = other.m_setorigin;
    m_setspan = other.m_setspan;
    m_setInfLimits = other.m_setInfLimits;
    m_setRefSys = other.m_setRefSys;
    m_isBuild = other.m_isBuild;
 
};
 
UStructMesh & UStructMesh::operator=(UStructMesh other){
    swap(other);
    return  *this;
};
 
void UStructMesh::swap(UStructMesh & x) noexcept
{
        // Number of cells
        std::swap(m_nx, x.m_nx);
        std::swap(m_ny, x.m_ny);
        std::swap(m_nz, x.m_nz);
 
        // Cell Spacing
        std::swap(m_dx, x.m_dx);
        std::swap(m_dy, x.m_dy);
        std::swap(m_dz, x.m_dz);
 
        // Copy cell edges and cell centers ----------------------------------------- //
        std::swap(m_xnode, x.m_xnode);
        std::swap(m_ynode, x.m_ynode);
        std::swap(m_znode, x.m_znode);
        std::swap(m_xedge, x.m_xedge);
        std::swap(m_yedge, x.m_yedge);
        std::swap(m_zedge, x.m_zedge);
 
        std::swap(m_origin_temp,x.m_origin_temp);
        std::swap(m_span_temp,x.m_span_temp);
        std::swap(m_inflimits_temp,x.m_inflimits_temp);
        std::swap(m_refsystem_temp,x.m_refsystem_temp);
        std::swap(m_shapetype_temp,x.m_shapetype_temp);
 
        std::swap(m_setorigin,x.m_setorigin);
        std::swap(m_setspan,x.m_setspan);
        std::swap(m_setInfLimits,x.m_setInfLimits);
        std::swap(m_setRefSys,x.m_setRefSys);
        std::swap(m_isBuild,x.m_isBuild);
 
        if(x.m_shape){
            std::unique_ptr<BasicShape> tempshape;
            switch(x.m_shape->getShapeType()){
                case ShapeType::WEDGE :
                {
                    const Wedge * pp = dynamic_cast<const Wedge*>(x.m_shape.get());
                    if (pp != nullptr)  tempshape = std::unique_ptr<BasicShape>(new Wedge(*(pp)));
                }
                break;
                case ShapeType::CYLINDER :
                {
                    const Cylinder * pp = dynamic_cast<const Cylinder*>(x.m_shape.get());
                    if (pp != nullptr)  tempshape = std::unique_ptr<BasicShape>(new Cylinder(*(pp)));
                }
                break;
                case ShapeType::SPHERE :
                {
                    const Sphere * pp = dynamic_cast<const Sphere*>(x.m_shape.get());
                    if (pp != nullptr)  tempshape = std::unique_ptr<BasicShape>(new Sphere(*(pp)));
                }
                break;
                default://CUBE
                {
                    const Cube * pp = dynamic_cast<const Cube*>(x.m_shape.get());
                    if (pp != nullptr)  tempshape = std::unique_ptr<BasicShape>(new Cube(*(pp)));
                }
                break;
            }
 
            x.m_shape = std::move(m_shape);
            m_shape = std::move(tempshape);
        }
}
 
 
const BasicShape * UStructMesh::getShape() const {
    return(m_shape.get());
}
 
BasicShape * UStructMesh::getShape(){
    return(m_shape.get());
}
 
darray3E UStructMesh::getOrigin(){
    if (getShape() == nullptr) return(m_origin_temp);
    return(getShape()->getOrigin());
}
 
darray3E UStructMesh::getSpan(){
    if (getShape() == nullptr) return(m_span_temp);
    return(getShape()->getSpan());
}
 
darray3E UStructMesh::getInfLimits(){
    if (getShape() == nullptr) return(m_inflimits_temp);
    return(getShape()->getInfLimits());
}
 
dmatrix33E UStructMesh::getRefSystem(){
    if (getShape() == nullptr) return(m_refsystem_temp);
    return(getShape()->getRefSystem());
}
 
darray3E UStructMesh::getScaling(){
    if (getShape() == nullptr) return(darray3E{{1,1,1}});
    return(getShape()->getScaling());
}
 
darray3E UStructMesh::getLocalSpan(){
    if (getShape() == nullptr) return(darray3E{{1,1,1}});
    return(getShape()->getLocalSpan());
}
 
ShapeType UStructMesh::getShapeType(){
    if (getShape() == nullptr) return(m_shapetype_temp);
    return(getShape()->getShapeType());
}
 
CoordType UStructMesh::getCoordTypex(){
    if (getShape() == nullptr) return(CoordType::CLAMPED);
    return(getShape()->getCoordinateType(0));
}
 
CoordType UStructMesh::getCoordType(int i){
    if (getShape() == nullptr) return(CoordType::CLAMPED);
    return(getShape()->getCoordinateType(i));
}
 
CoordType UStructMesh::getCoordTypey(){
    if (getShape() == nullptr) return(CoordType::CLAMPED);
    return(getShape()->getCoordinateType(1));
}
 
CoordType UStructMesh::getCoordTypez(){
    if (getShape() == nullptr) return(CoordType::CLAMPED);
    return(getShape()->getCoordinateType(2));
}
 
std::array<CoordType, 3> UStructMesh::getCoordType(){
    std::array<CoordType, 3> types;
    for (int i=0; i<3; i++){
        types[i] = getShape()->getCoordinateType(i);
    }
    return(types);
}
 
darray3E UStructMesh::getSpacing(){
    darray3E res;
    darray3E scale = getScaling();
    res[0] = m_dx*scale[0]; res[1] =m_dy*scale[1]; res[2] = m_dz*scale[2];
    return(res);
};
 
iarray3E UStructMesh::getDimension(){
 
    iarray3E res;
    res[0] = m_nx + 1;
    res[1] = m_ny + 1;
    res[2] = m_nz + 1;
    return(res);
};
 
darray3E UStructMesh::getLocalCCell(int index){
 
    int i,j,k;
    accessCellIndex(index, i,j,k);
    darray3E res = getLocalCCell(i,j,k);
 
    return(res);
};
 
darray3E UStructMesh::getLocalCCell(int i_, int j_, int k_){
 
    darray3E res;
    res[0] = m_xnode[i_];
    res[1] = m_ynode[j_];
    res[2] = m_znode[k_];
 
    return(res);
};
 
darray3E UStructMesh::getLocalPoint(int index){
 
    int i,j,k;
    accessPointIndex(index, i,j,k);
    darray3E res = getLocalPoint(i,j,k);
 
    return(res);
};
darray3E UStructMesh::getLocalPoint(int i_, int j_, int k_){
 
    darray3E res;
    res[0] = m_xedge[i_];
    res[1] = m_yedge[j_];
    res[2] = m_zedge[k_];
 
    return(res);
};
 
darray3E UStructMesh::getGlobalCCell(int index){
 
    darray3E res = getLocalCCell(index);
    return(transfToGlobal(res));
};
darray3E UStructMesh::getGlobalCCell(int i_, int j_, int k_){
 
    darray3E res = getLocalCCell(i_,j_,k_);
    return(transfToGlobal(res));
};
 
darray3E UStructMesh::getGlobalPoint(int index){
    darray3E res = getLocalPoint(index);
    return(transfToGlobal(res));
};
darray3E UStructMesh::getGlobalPoint(int i_, int j_, int k_){
 
    darray3E res = getLocalPoint(i_,j_,k_);
    return(transfToGlobal(res));
};
 
ivector1D UStructMesh::getCellNeighs(int i_, int j_, int k_){
 
    ivector1D result(8);
    result[0] = accessPointIndex(i_, j_, k_);
    result[1] = accessPointIndex(i_+1, j_, k_);
    result[2] = accessPointIndex(i_+1, j_+1, k_);
    result[3] = accessPointIndex(i_, j_+1, k_);
    result[4] = accessPointIndex(i_, j_, k_+1);
    result[5] = accessPointIndex(i_+1, j_, k_+1);
    result[6] = accessPointIndex(i_+1, j_+1, k_+1);
    result[7] = accessPointIndex(i_, j_+1, k_+1);
    return(result);
}
 
ivector1D UStructMesh::getCellNeighs(int index){
 
    ivector1D pp(3,0);
    accessCellIndex(index, pp[0],pp[1], pp[2]);
    return(getCellNeighs(pp[0],pp[1],pp[2]));
}
 
dvecarr3E UStructMesh::getLocalCoords(){
 
    int np = (m_nx+1)*(m_ny+1)*(m_nz+1);
    dvecarr3E coords(np);
    for (int i=0; i<np; i++){
        coords[i] = getLocalPoint(i);
    }
    return coords;
};
 
dvecarr3E UStructMesh::getGlobalCoords(){
    int np = (m_nx+1)*(m_ny+1)*(m_nz+1);
    dvecarr3E coords(np);
    for (int i=0; i<np; i++){
        coords[i] = getGlobalPoint(i);
    }
    return coords;
};
 
dvecarr3E UStructMesh::getLocalCellCentroids(){
 
    int np = (m_nx)*(m_ny)*(m_nz);
    dvecarr3E coords(np);
    for (int i=0; i<np; i++){
        coords[i] = getLocalCCell(i);
    }
    return coords;
};
 
dvecarr3E UStructMesh::getGlobalCellCentroids(){
    int np = (m_nx)*(m_ny)*(m_nz);
    dvecarr3E coords(np);
    for (int i=0; i<np; i++){
        coords[i] = getGlobalCCell(i);
    }
    return coords;
};
 
void UStructMesh::setOrigin(darray3E origin){
    if (getShape() == nullptr){
        m_origin_temp = origin;
        m_setorigin = true;
    }else{
        getShape()->setOrigin(origin);
    }
}
 
void UStructMesh::setSpan(double s0, double s1, double s2){
 
    if (getShape() == nullptr){
        m_span_temp[0] = s0;
        m_span_temp[1] = s1;
        m_span_temp[2] = s2;
        m_setspan = true;
    }else{
        getShape()->setSpan( s0, s1,s2);
    }
    m_isBuild = false;
}
 
void UStructMesh::setSpan(darray3E s){
    setSpan( s[0], s[1], s[2]);
}
 
void UStructMesh::setInfLimits(double inflim, int dir){
    if (getShape() == nullptr){
        m_inflimits_temp[dir] = inflim;
        m_setInfLimits = true;
    }else{
        getShape()->setInfLimits(inflim, dir);
    }
}
 
void UStructMesh::setInfLimits(darray3E inflim){
    setInfLimits(inflim[0], 0);
    setInfLimits(inflim[1], 1);
    setInfLimits(inflim[2], 2);
}
 
void UStructMesh::setRefSystem(darray3E axis0, darray3E axis1, darray3E axis2){
 
    if (getShape() == nullptr){
        m_refsystem_temp[0] = axis0;
        m_refsystem_temp[1] = axis1;
        m_refsystem_temp[2] = axis2;
        m_setRefSys = true;
    }else{
        getShape()->setRefSystem(axis0, axis1, axis2);
    }
}
 
void UStructMesh::setRefSystem(int label, darray3E axis){
 
    if(getShape() == nullptr){
        BasicShape * temp = new Cube();
        temp->setRefSystem(label, axis);
        m_refsystem_temp = temp->getRefSystem();
        delete temp;
        temp = nullptr;
        m_setRefSys = true;
    }else{
        getShape()->setRefSystem(label,axis);
    }
}
 
void UStructMesh::setRefSystem(dmatrix33E axes){
    if (getShape() == nullptr){
 
        m_refsystem_temp[0] = axes[0];
        m_refsystem_temp[1] = axes[1];
        m_refsystem_temp[2] = axes[2];
        m_setRefSys = true;
    }else{
        getShape()->setRefSystem(axes);
    }
}
 
 void UStructMesh::setDimension(ivector1D dim){
     m_nx = dim[0] - 1;
     m_ny = dim[1] - 1;
     m_nz = dim[2] - 1;
 
     m_isBuild = false;
};
 
  void UStructMesh::setDimension(iarray3E dim){
     m_nx = dim[0] - 1;
     m_ny = dim[1] - 1;
     m_nz = dim[2] - 1;
 
     m_isBuild = false;
 };
 
  void UStructMesh::setShape(int itype){
      ShapeType type = static_cast<ShapeType>(itype);
      UStructMesh::setShape(type);
  }
 
  void UStructMesh::setShape(ShapeType type){
      //create internal shape using unique_ptr member.
      // unlink external shape eventually
        if(m_shape){
 
            m_origin_temp = m_shape.get()->getOrigin();
            m_span_temp = m_shape.get()->getSpan();
            m_inflimits_temp = m_shape.get()->getInfLimits();
            m_refsystem_temp = m_shape.get()->getRefSystem();
            m_setorigin = true;
            m_setspan = true;
            m_setInfLimits = true;
            m_setRefSys = true;
 
            m_shape.reset(nullptr);
        }
 
        darray3E origin{{0,0,0}};
        dmatrix33E spanMat;
        spanMat[0].fill(1); spanMat[2].fill(1); spanMat[2].fill(1);
        spanMat[1][1] = spanMat[2][1] = 2*BITPIT_PI;
        spanMat[2][2] = BITPIT_PI;
 
        if(m_setorigin){origin = m_origin_temp; m_setorigin= false;}
        if(m_setspan){
            for(int i=0; i<3; ++i) spanMat[i] = m_span_temp;
            m_setspan = false;
        }
 
        switch(type){
            case ShapeType::WEDGE :
                m_shape = std::unique_ptr<BasicShape>(new Wedge(origin,spanMat[0]));
                break;
            case ShapeType::CYLINDER :
                m_shape = std::unique_ptr<BasicShape>(new Cylinder(origin,spanMat[1]));
                break;
            case ShapeType::SPHERE :
                m_shape = std::unique_ptr<BasicShape>(new Sphere(origin, spanMat[2]));
                break;
            default://CUBE
                m_shape = std::unique_ptr<BasicShape>(new Cube(origin, spanMat[0]));
                break;
        }
 
        if(m_setInfLimits){
            m_shape.get()->setInfLimits(m_inflimits_temp[0],0);
            m_shape.get()->setInfLimits(m_inflimits_temp[1],1);
            m_shape.get()->setInfLimits(m_inflimits_temp[2],2);
            m_setInfLimits = false;
        }
 
        if(m_setRefSys){
            m_shape.get()->setRefSystem(m_refsystem_temp);
            m_setRefSys = false;
        }
 
        m_isBuild = false;
}
 
void UStructMesh::setShape(const BasicShape * shape){
 
    if(shape == nullptr) return;
    m_shape.reset(nullptr);
    m_setorigin = false;
    m_setspan = false;
    m_setInfLimits = false;
    m_setRefSys = false;
 
    switch(shape->getShapeType()){
        case ShapeType::WEDGE :
        {
            const Wedge * pp = dynamic_cast<const Wedge*>(shape);
            if (pp != nullptr)  m_shape = std::unique_ptr<BasicShape>(new Wedge(*(pp)));
        }
        break;
        case ShapeType::CYLINDER :
            {
                const Cylinder * pp = dynamic_cast<const Cylinder*>(shape);
                if (pp != nullptr)  m_shape = std::unique_ptr<BasicShape>(new Cylinder(*(pp)));
            }
            break;
        case ShapeType::SPHERE :
            {
                const Sphere * pp = dynamic_cast<const Sphere*>(shape);
                if (pp != nullptr)  m_shape = std::unique_ptr<BasicShape>(new Sphere(*(pp)));
            }
            break;
        default://CUBE
            {
                const Cube * pp = dynamic_cast<const Cube*>(shape);
                if (pp != nullptr)  m_shape = std::unique_ptr<BasicShape>(new Cube(*(pp)));
            }
            break;
    }
 
    m_isBuild = false;
};
 
 
void UStructMesh::setCoordTypex(CoordType type){
    if (getShape() == nullptr) return;
    getShape()->setCoordinateType(type,0);
}
 
void UStructMesh::setCoordType(CoordType type, int i){
    if (getShape() == nullptr) return;
    getShape()->setCoordinateType(type,i);
}
 
void UStructMesh::setCoordTypey(CoordType type){
    if (getShape() == nullptr) return;
    getShape()->setCoordinateType(type,1);
}
 
void UStructMesh::setCoordTypez(CoordType type){
    if (getShape() == nullptr) return;
    getShape()->setCoordinateType(type,2);
}
 
void UStructMesh::setCoordType(std::array<CoordType, 3> types){
    for (int i=0; i<3; i++){
        getShape()->setCoordinateType(types[i],i);
    }
}
 
void UStructMesh::setMesh(darray3E & origin, darray3E &span, ShapeType type, iarray3E & dimensions){
 
    if(m_shape){m_shape.reset(nullptr);}
 
    setShape(type);
    setOrigin(origin);
    setSpan(span);
    setDimension(dimensions);
    build();
};
 
 
void UStructMesh::setMesh(darray3E & origin, darray3E &span, ShapeType type, dvector1D & spacing){
 
    ivector1D dimLimit(3,2);
    //create internal shape using unique_ptr member.
    if(m_shape){m_shape.reset(nullptr);}
 
    switch(type){
        case ShapeType::CYLINDER :
            dimLimit[1] = 4;
            break;
        case ShapeType::SPHERE :
            dimLimit[1] = 5; dimLimit[2] = 3;
            break;
        default://CUBE //WEDGE
            break;
    }
 
    setShape(type);
    setOrigin(origin);
    setSpan(span);
 
    darray3E span2 = getSpan();
    iarray3E dim;
 
    for(int i=0; i<3; ++i){
        if(spacing[i] != 0.0) {
            dim[i] = (int) std::floor(span2[i]/spacing[i] +0.5) + 1;
        }else{
            dim[i] = dimLimit[i];
        }
    }
 
    setDimension(dim);
    build();
};
 
void UStructMesh::setMesh(BasicShape * shape, iarray3E & dimensions){
 
    if(m_shape){m_shape.reset(nullptr);}
 
    setShape(shape);
    setDimension(dimensions);
    build();
};
 
void UStructMesh::setMesh(BasicShape * shape, dvector1D & spacing){
 
    ivector1D dimLimit(3,2);
    //create internal shape using unique_ptr member.
    if(m_shape){m_shape.reset(nullptr);}
 
    switch(shape->getShapeType()){
        case ShapeType::CYLINDER :
            dimLimit[1] = 4;
            break;
        case ShapeType::SPHERE :
            dimLimit[1] = 5; dimLimit[2] = 3;
            break;
        default://CUBE //WEDGE
            break;
    }
 
    setShape(shape);
 
    darray3E span2 = getSpan();
    iarray3E dim;
 
    for(int i=0; i<3; ++i){
        if(spacing[i] != 0.0) {
            dim[i] = (int) std::floor(span2[i]/spacing[i] +0.5) + 1;
        }else{
            dim[i] = dimLimit[i];
        }
    }
 
    setDimension(dim);
    build();
};
 
 
void UStructMesh::clearMesh(){
 
    m_shape.reset(nullptr);
    m_nx=0; m_ny=0; m_nz=0;
    m_dx=0.0; m_dy=0.0; m_dz=0.0;
 
 
    m_setorigin = false;
    m_setspan = false;
    m_setInfLimits = false;
    m_setRefSys = false;
    m_isBuild = false;
 
    m_origin_temp = {{0.0,0.0,0.0}};
    m_span_temp = {{1.0,1.0,1.0}};
    m_inflimits_temp = {{0.0,0.0,0.0}};
    for(int i=0; i<3; ++i){m_refsystem_temp[i].fill(0.0); m_refsystem_temp[i][i] = 1.0;}
    m_shapetype_temp = ShapeType::CUBE;
 
    destroyNodalStructure();
};
 
 
void UStructMesh::locateCellByPoint(darray3E & point, int &i, int &j, int &k){
 
    darray3E P = transfToLocal(point);
    darray3E locOr = getShape()->getLocalOrigin();
 
    i = std::min(m_nx-1, std::max(0, (int) std::floor((P[0]-locOr[0])/m_dx)));
    j = std::min(m_ny-1, std::max(0, (int) std::floor((P[1]-locOr[1])/m_dy)));
    k = std::min(m_nz-1, std::max(0, (int) std::floor((P[2]-locOr[2])/m_dz)));
};
 
void UStructMesh::locateCellByPoint(dvector1D & point, int &i, int &j, int &k){
        darray3E temp = conArray<double,3>(point);
        locateCellByPoint(temp,i,j,k);
};
 
int UStructMesh::accessCellIndex(int i, int j, int k){
 
    int index = m_ny * m_nz * i + m_nz * j + k;
    return(index);
};
 
void UStructMesh::accessCellIndex(int N_, int & i, int & j, int & k){
    k = N_ % m_nz;
    int index = N_ / m_nz;
    j = index % m_ny;
    i = index / m_ny;
};
 
void UStructMesh::accessPointIndex(int N_, int &i, int &j, int &k){
 
    k = N_ % (m_nz+1);
    int index = N_ / (m_nz+1);
    j = index % (m_ny+1);
    i = index / (m_ny+1);
};
 
darray3E    UStructMesh::transfToGlobal( darray3E & point){
    return(getShape()->toWorldCoord(point));
};
 
dvector1D   UStructMesh::transfToGlobal( dvector1D & point){
    darray3E temp = conArray<double,3>(point);
    darray3E temp2 = getShape()->toWorldCoord(temp);
    return(conVect(temp2));
};
 
dvecarr3E   UStructMesh::transfToGlobal( dvecarr3E & list_points){
    int size = list_points.size();
    dvecarr3E result(size);
    for(int i=0; i<size; ++i){
        result[i] = transfToGlobal(list_points[i]);
    }
    return(result);
};
 
darray3E    UStructMesh::transfToLocal( darray3E & point){
    return(getShape()->toLocalCoord(point));
};
 
dvector1D   UStructMesh::transfToLocal( dvector1D & point){
    darray3E temp = conArray<double,3>(point);
    darray3E temp2 = getShape()->toLocalCoord(temp);
    return(conVect(temp2));
 
};
 
dvecarr3E   UStructMesh::transfToLocal( dvecarr3E & list_points){
    int size = list_points.size();
    dvecarr3E result(size);
    for(int i=0; i<size; ++i){
        result[i] = transfToLocal(list_points[i]);
    }
    return(result);
 
};
 
double UStructMesh::interpolateCellData(darray3E & point, dvector1D & celldata){
 
    int i0, j0, k0, ip, jp, kp;
    double wx0,wx1,wy0,wy1,wz0,wz1;
 
    locateCellByPoint(point, i0, j0, k0);
    darray3E P = transfToLocal(point);
    if (P[0] > m_xnode[i0])     { ip = std::min(i0+1, m_nx-1);   }
    else                    { ip = std::max(0, i0-1);      }
    if (P[1] > m_ynode[j0])     { jp = std::min(j0+1, m_ny-1);   }
    else                    { jp = std::max(0, j0-1);      }
    if (P[2] > m_znode[k0])     { kp = std::min(k0+1, m_nz-1);   }
    else                    { kp = std::max(0, k0-1);      }
 
    // Interpolation weights
    wx1 = std::max(0.0, std::min(1.0, std::abs((P[0] - m_xnode[i0])/m_dx)));     wx0 = 1.0 - wx1;
    wy1 = std::max(0.0, std::min(1.0, std::abs((P[1] - m_ynode[j0])/m_dy)));     wy0 = 1.0 - wy1;
    wz1 = std::max(0.0, std::min(1.0, std::abs((P[2] - m_znode[k0])/m_dz)));     wz0 = 1.0 - wz1;
 
    // Interpolation
    double result  = wz0 * wx0 * wy0 * celldata[accessCellIndex(i0,j0,k0)]
    + wz0 * wx0 * wy1 * celldata[accessCellIndex(i0,jp,k0)]
    + wz0 * wx1 * wy0 * celldata[accessCellIndex(ip,j0,k0)]
    + wz0 * wx1 * wy1 * celldata[accessCellIndex(ip,jp,k0)]
    + wz1 * wx0 * wy0 * celldata[accessCellIndex(i0,j0,kp)]
    + wz1 * wx0 * wy1 * celldata[accessCellIndex(i0,jp,kp)]
    + wz1 * wx1 * wy0 * celldata[accessCellIndex(ip,j0,kp)]
    + wz1 * wx1 * wy1 * celldata[accessCellIndex(ip,jp,kp)];
 
    return(result);
};
 
int UStructMesh::interpolateCellData(darray3E & point, ivector1D & celldata){
 
    int i0, j0, k0, ip, jp, kp;
    double wx0,wx1,wy0,wy1,wz0,wz1;
 
    locateCellByPoint(point, i0, j0, k0);
    darray3E P = transfToLocal(point);
    if (P[0] > m_xnode[i0])     { ip = std::min(i0+1, m_nx-1);   }
    else                    { ip = std::max(0, i0-1);      }
    if (P[1] > m_ynode[j0])     { jp = std::min(j0+1, m_ny-1);   }
    else                    { jp = std::max(0, j0-1);      }
    if (P[2] > m_znode[k0])     { kp = std::min(k0+1, m_nz-1);   }
    else                    { kp = std::max(0, k0-1);      }
 
    // Interpolation weights
    wx1 = std::max(0.0, std::min(1.0, std::abs((P[0] - m_xnode[i0])/m_dx)));     wx0 = 1.0 - wx1;
    wy1 = std::max(0.0, std::min(1.0, std::abs((P[1] - m_ynode[j0])/m_dy)));     wy0 = 1.0 - wy1;
    wz1 = std::max(0.0, std::min(1.0, std::abs((P[2] - m_znode[k0])/m_dz)));     wz0 = 1.0 - wz1;
 
    // Interpolation
    double result  =    wz0 * wx0 * wy0 * (double)celldata[accessCellIndex(i0,j0,k0)]
    + wz0 * wx0 * wy1 * (double)celldata[accessCellIndex(i0,jp,k0)]
    + wz0 * wx1 * wy0 *(double) celldata[accessCellIndex(ip,j0,k0)]
    + wz0 * wx1 * wy1 *(double) celldata[accessCellIndex(ip,jp,k0)]
    + wz1 * wx0 * wy0 *(double) celldata[accessCellIndex(i0,j0,kp)]
    + wz1 * wx0 * wy1 *(double) celldata[accessCellIndex(i0,jp,kp)]
    + wz1 * wx1 * wy0 *(double) celldata[accessCellIndex(ip,j0,kp)]
    + wz1 * wx1 * wy1 *(double) celldata[accessCellIndex(ip,jp,kp)];
 
    int result2 =std::floor(result+0.5);
    return(result2);
};
 
darray3E UStructMesh::interpolateCellData(darray3E & point, dvecarr3E & celldata){
 
    int i0, j0, k0, ip, jp, kp;
    double wx0,wx1,wy0,wy1,wz0,wz1;
 
    locateCellByPoint(point, i0, j0, k0);
    darray3E P = transfToLocal(point);
 
    if (P[0] > m_xnode[i0])     { ip = std::min(i0+1, m_nx-1);   }
    else                    { ip = std::max(0, i0-1);      }
    if (P[1] > m_ynode[j0])     { jp = std::min(j0+1, m_ny-1);   }
    else                    { jp = std::max(0, j0-1);      }
    if (P[2] > m_znode[k0])     { kp = std::min(k0+1, m_nz-1);   }
    else                    { kp = std::max(0, k0-1);      }
 
    // Interpolation weights
    wx1 = std::max(0.0, std::min(1.0, std::abs((P[0] - m_xnode[i0])/m_dx)));     wx0 = 1.0 - wx1;
    wy1 = std::max(0.0, std::min(1.0, std::abs((P[1] - m_ynode[j0])/m_dy)));     wy0 = 1.0 - wy1;
    wz1 = std::max(0.0, std::min(1.0, std::abs((P[2] - m_znode[k0])/m_dz)));     wz0 = 1.0 - wz1;
 
    // Interpolation
    darray3E result = wz0 * wx0 * wy0 * celldata[accessCellIndex(i0,j0,k0)]
    + wz0 * wx0 * wy1 * celldata[accessCellIndex(i0,jp,k0)]
    + wz0 * wx1 * wy0 * celldata[accessCellIndex(ip,j0,k0)]
    + wz0 * wx1 * wy1 * celldata[accessCellIndex(ip,jp,k0)]
    + wz1 * wx0 * wy0 * celldata[accessCellIndex(i0,j0,kp)]
    + wz1 * wx0 * wy1 * celldata[accessCellIndex(i0,jp,kp)]
    + wz1 * wx1 * wy0 * celldata[accessCellIndex(ip,j0,kp)]
    + wz1 * wx1 * wy1 * celldata[accessCellIndex(ip,jp,kp)];
 
    return(result);
};
 
double UStructMesh::interpolatePointData(darray3E & point, dvector1D & pointdata){
 
    int i0, j0, k0, ip, jp, kp;
    double wx0,wx1,wy0,wy1,wz0,wz1;
 
    darray3E P = transfToLocal(point);
    darray3E locOr = getShape()->getLocalOrigin();
    i0 = std::max(0, std::min(m_nx, (int) std::floor((P[0]-locOr[0])/m_dx)));
    j0 = std::max(0, std::min(m_ny, (int) std::floor((P[1]-locOr[1])/m_dy)));
    k0 = std::max(0, std::min(m_nz, (int) std::floor((P[2]-locOr[2])/m_dz)));
 
    if (P[0] >= m_xedge[i0])    { ip = std::min(i0+1, m_nx);   }
    else                    { ip = std::max(0, i0-1);      }
    if (P[1] >= m_yedge[j0])    { jp = std::min(j0+1, m_ny);   }
    else                    { jp = std::max(0, j0-1);      }
    if (P[2] >= m_zedge[k0])    { kp = std::min(k0+1, m_nz);   }
    else                    { kp = std::max(0, k0-1);      }
 
    // Interpolation weights
    wx1 = std::max(0.0, std::min(1.0, std::abs((P[0] - m_xedge[i0])/m_dx)));     wx0 = 1.0 - wx1;
    wy1 = std::max(0.0, std::min(1.0, std::abs((P[1] - m_yedge[j0])/m_dy)));     wy0 = 1.0 - wy1;
    wz1 = std::max(0.0, std::min(1.0, std::abs((P[2] - m_zedge[k0])/m_dz)));     wz0 = 1.0 - wz1;
 
    // Interpolation
    double result =  wz0 * wx0 * wy0 * pointdata[accessPointIndex(i0,j0,k0)]
    + wz0 * wx0 * wy1 * pointdata[accessPointIndex(i0,jp,k0)]
    + wz0 * wx1 * wy0 * pointdata[accessPointIndex(ip,j0,k0)]
    + wz0 * wx1 * wy1 * pointdata[accessPointIndex(ip,jp,k0)]
    + wz1 * wx0 * wy0 * pointdata[accessPointIndex(i0,j0,kp)]
    + wz1 * wx0 * wy1 * pointdata[accessPointIndex(i0,jp,kp)]
    + wz1 * wx1 * wy0 * pointdata[accessPointIndex(ip,j0,kp)]
    + wz1 * wx1 * wy1 * pointdata[accessPointIndex(ip,jp,kp)];
    return(result);
};
 
int UStructMesh::interpolatePointData(darray3E & point, ivector1D & pointdata){
 
    int i0, j0, k0, ip, jp, kp;
    double wx0,wx1,wy0,wy1,wz0,wz1;
 
    darray3E P = transfToLocal(point);
    darray3E locOr = getShape()->getLocalOrigin();
    i0 = std::max(0, std::min(m_nx, (int) std::floor((P[0]-locOr[0])/m_dx)));
    j0 = std::max(0, std::min(m_ny, (int) std::floor((P[1]-locOr[1])/m_dy)));
    k0 = std::max(0, std::min(m_nz, (int) std::floor((P[2]-locOr[2])/m_dz)));
 
    if (P[0] >= m_xedge[i0])    { ip = std::min(i0+1, m_nx);   }
    else                    { ip = std::max(0, i0-1);      }
    if (P[1] >= m_yedge[j0])    { jp = std::min(j0+1, m_ny);   }
    else                    { jp = std::max(0, j0-1);      }
    if (P[2] >= m_zedge[k0])    { kp = std::min(k0+1, m_nz);   }
    else                    { kp = std::max(0, k0-1);      }
 
    // Interpolation weights
    wx1 = std::max(0.0, std::min(1.0, std::abs((P[0] - m_xedge[i0])/m_dx)));     wx0 = 1.0 - wx1;
    wy1 = std::max(0.0, std::min(1.0, std::abs((P[1] - m_yedge[j0])/m_dy)));     wy0 = 1.0 - wy1;
    wz1 = std::max(0.0, std::min(1.0, std::abs((P[2] - m_zedge[k0])/m_dz)));     wz0 = 1.0 - wz1;
 
    // Interpolation
    double result = wz0 * wx0 * wy0 * (double)pointdata[accessPointIndex(i0,j0,k0)]
    + wz0 * wx0 * wy1 * (double)pointdata[accessPointIndex(i0,jp,k0)]
    + wz0 * wx1 * wy0 *(double) pointdata[accessPointIndex(ip,j0,k0)]
    + wz0 * wx1 * wy1 *(double) pointdata[accessPointIndex(ip,jp,k0)]
    + wz1 * wx0 * wy0 *(double) pointdata[accessPointIndex(i0,j0,kp)]
    + wz1 * wx0 * wy1 *(double) pointdata[accessPointIndex(i0,jp,kp)]
    + wz1 * wx1 * wy0 *(double) pointdata[accessPointIndex(ip,j0,kp)]
    + wz1 * wx1 * wy1 *(double) pointdata[accessPointIndex(ip,jp,kp)];
 
    int result2 = std::floor(result+0.5);
    return(result2);
};
 
darray3E UStructMesh::interpolatePointData(darray3E & point, dvecarr3E & pointdata){
 
    int i0, j0, k0, ip, jp, kp;
    double wx0,wx1,wy0,wy1,wz0,wz1;
 
    darray3E P = transfToLocal(point);
    darray3E locOr = getShape()->getLocalOrigin();
    i0 = std::max(0, std::min(m_nx, (int) std::floor((P[0]-locOr[0])/m_dx)));
    j0 = std::max(0, std::min(m_ny, (int) std::floor((P[1]-locOr[1])/m_dy)));
    k0 = std::max(0, std::min(m_nz, (int) std::floor((P[2]-locOr[2])/m_dz)));
 
 
    if (P[0] >= m_xedge[i0])    { ip = std::min(i0+1, m_nx);   }
    else                    { ip = std::max(0, i0-1);      }
    if (P[1] >= m_yedge[j0])    { jp = std::min(j0+1, m_ny);   }
    else                    { jp = std::max(0, j0-1);      }
    if (P[2] >= m_zedge[k0])    { kp = std::min(k0+1, m_nz);   }
    else                    { kp = std::max(0, k0-1);      }
 
    // Interpolation weights
    wx1 = std::max(0.0, std::min(1.0, std::abs((P[0] - m_xedge[i0])/m_dx)));     wx0 = 1.0 - wx1;
    wy1 = std::max(0.0, std::min(1.0, std::abs((P[1] - m_yedge[j0])/m_dy)));     wy0 = 1.0 - wy1;
    wz1 = std::max(0.0, std::min(1.0, std::abs((P[2] - m_zedge[k0])/m_dz)));     wz0 = 1.0 - wz1;
 
    // Interpolation
    darray3E result = wz0 * wx0 * wy0 * pointdata[accessPointIndex(i0,j0,k0)]
    + wz0 * wx0 * wy1 * pointdata[accessPointIndex(i0,jp,k0)]
    + wz0 * wx1 * wy0 * pointdata[accessPointIndex(ip,j0,k0)]
    + wz0 * wx1 * wy1 * pointdata[accessPointIndex(ip,jp,k0)]
    + wz1 * wx0 * wy0 * pointdata[accessPointIndex(i0,j0,kp)]
    + wz1 * wx0 * wy1 * pointdata[accessPointIndex(i0,jp,kp)]
    + wz1 * wx1 * wy0 * pointdata[accessPointIndex(ip,j0,kp)]
    + wz1 * wx1 * wy1 * pointdata[accessPointIndex(ip,jp,kp)];
 
    return(result);
};
 
 
void UStructMesh::plotCloud( std::string & folder , std::string outfile,
                             int counterFile, bool codexFlag,
                             const ivector1D & labels,
                             dvecarr3E * extPoints){
 
    bitpit::VTKFormat codex = bitpit::VTKFormat::ASCII;
    if(codexFlag){codex=bitpit::VTKFormat::APPENDED;}
 
    iarray3E dim = getDimension();
    int sizeTot = dim[0]*dim[1]*dim[2];
 
    dvecarr3E activeP;
    if(extPoints != nullptr && (int)extPoints->size() == sizeTot){
        activeP = *extPoints;
    }else{
        activeP.resize(sizeTot);
        for(int i=0; i<sizeTot; i++){
            activeP[i] = getGlobalPoint(i);
        }
    }
 
    ivector1D conn(activeP.size());
    {
        int counter = 0;
        for(auto & val: conn){
            val = counter;
            ++counter;
        }
    }
    ivector1D datalab = labels;
    datalab.resize(activeP.size(), -1);
 
    bitpit::VTKUnstructuredGrid vtk(folder, outfile, bitpit::VTKElementType::VERTEX);
    vtk.setGeomData( bitpit::VTKUnstructuredField::POINTS, activeP) ;
    vtk.setGeomData( bitpit::VTKUnstructuredField::CONNECTIVITY, conn) ;
    vtk.setDimensions(conn.size(), activeP.size());
    vtk.addData("labels",bitpit::VTKFieldType::SCALAR, bitpit::VTKLocation::POINT, datalab) ;
    vtk.setCodex(codex);
    if(counterFile>=0){vtk.setCounter(counterFile);}
 
    vtk.write();
 
};
 
void UStructMesh::plotCloudScalar( std::string folder , std::string outfile, int counterFile, bool codexFlag, dvector1D & data){
 
    bitpit::VTKFormat codex = bitpit::VTKFormat::ASCII;
    if(codexFlag){codex=bitpit::VTKFormat::APPENDED;}
 
    iarray3E dim = getDimension();
    int sizeTot = dim[0]*dim[1]*dim[2];
 
    dvecarr3E activeP(sizeTot);
    for(int i=0; i<sizeTot; i++){
        activeP[i] = getGlobalPoint(i);
    }
 
    ivector1D conn(activeP.size());
    {
        int counter = 0;
        for(auto & val: conn){
            val = counter;
            ++counter;
        }
    }
 
    bitpit::VTKUnstructuredGrid vtk(folder, outfile, bitpit::VTKElementType::VERTEX);
    vtk.setGeomData( bitpit::VTKUnstructuredField::POINTS, activeP) ;
    vtk.setGeomData( bitpit::VTKUnstructuredField::CONNECTIVITY, conn) ;
    vtk.setDimensions(conn.size(), activeP.size());
    vtk.setCodex(codex);
    if(counterFile>=0){vtk.setCounter(counterFile);}
 
    bitpit::VTKLocation loc = bitpit::VTKLocation::POINT;
    vtk.addData("field",bitpit::VTKFieldType::SCALAR, loc, data) ;
    vtk.write();
};
 
void UStructMesh::plotGrid(std::string & folder, std::string outfile ,
                           int counterFile, bool codexFlag,
                           const ivector1D & labels,
                           dvecarr3E * extPoints)
{
    bitpit::VTKFormat codex = bitpit::VTKFormat::ASCII;
    if(codexFlag){codex=bitpit::VTKFormat::APPENDED;}
 
    iarray3E dim = getDimension();
    int sizePt = dim[0]*dim[1]*dim[2];
    int sizeCl = (dim[0]-1)*(dim[1]-1)*(dim[2]-1);
 
    dvecarr3E activeP(sizePt);
    ivector2D activeConn(sizeCl, ivector1D(8,0));
 
    if(extPoints != nullptr && (int)extPoints->size() == sizePt){
        activeP = *extPoints;
    }
    else{
        for(int i=0; i<sizePt; i++){
            activeP[i] = getGlobalPoint(i);
        }
    }
 
    for(int i=0; i<sizeCl; ++i){
        activeConn[i] = getCellNeighs(i);
    }
 
    ivector1D datalab = labels;
    datalab.resize(activeP.size(),-1);
 
    bitpit::VTKElementType elDM = bitpit::VTKElementType::HEXAHEDRON;
    bitpit::VTKUnstructuredGrid vtk(folder, outfile, elDM);
    vtk.setGeomData( bitpit::VTKUnstructuredField::POINTS, activeP) ;
    vtk.setGeomData( bitpit::VTKUnstructuredField::CONNECTIVITY, activeConn) ;
    vtk.setDimensions(sizeCl, sizePt);
    vtk.addData("labels",bitpit::VTKFieldType::SCALAR, bitpit::VTKLocation::POINT, datalab) ;
    vtk.setCodex(codex);
    if(counterFile>=0){vtk.setCounter(counterFile);}
 
    vtk.write();
 
};
 
void UStructMesh::plotGridScalar(std::string folder, std::string outfile , int counterFile, bool codexFlag, dvector1D & data){
 
    bitpit::VTKFormat codex = bitpit::VTKFormat::ASCII;
    if(codexFlag){codex=bitpit::VTKFormat::APPENDED;}
 
    iarray3E dim = getDimension();
    int sizePt = dim[0]*dim[1]*dim[2];
    int sizeCl = (dim[0]-1)*(dim[1]-1)*(dim[2]-1);
 
    dvecarr3E activeP(sizePt);
    ivector2D activeConn(sizeCl, ivector1D(8,0));
 
    for(int i=0; i<sizePt; i++){
        activeP[i] = getGlobalPoint(i);
    }
 
    for(int i=0; i<sizeCl; ++i){
        activeConn[i] = getCellNeighs(i);
    }
 
    bitpit::VTKElementType  elDM = bitpit::VTKElementType::HEXAHEDRON;
    bitpit::VTKUnstructuredGrid vtk(folder, outfile, elDM);
    vtk.setGeomData( bitpit::VTKUnstructuredField::POINTS, activeP) ;
    vtk.setGeomData( bitpit::VTKUnstructuredField::CONNECTIVITY, activeConn) ;
    vtk.setDimensions(sizeCl, sizePt);
 
    vtk.setCodex(codex);
    if(counterFile>=0){vtk.setCounter(counterFile);}
 
    bitpit::VTKLocation loc = bitpit::VTKLocation::POINT;
    if ((int)data.size() == sizeCl) loc = bitpit::VTKLocation::CELL;
 
    vtk.addData("field",bitpit::VTKFieldType::SCALAR, loc, data) ;
    vtk.write();
 
};
 
 
void UStructMesh::destroyNodalStructure(){
    m_xnode.clear();
    m_ynode.clear();
    m_znode.clear();
    m_xedge.clear();
    m_yedge.clear();
    m_zedge.clear();
};
 
void UStructMesh::reshapeNodalStructure(){
    destroyNodalStructure();
    resizeMesh();
};
 
void UStructMesh::resizeMesh(){
    // Cell centers
    m_xnode.resize(m_nx, 0.0);
    m_ynode.resize(m_ny, 0.0);
    m_znode.resize(m_nz, 0.0);
 
    // Points
    m_xedge.resize(m_nx+1, 0.0);
    m_yedge.resize(m_ny+1, 0.0);
    m_zedge.resize(m_nz+1, 0.0);
};
 
bool UStructMesh::build(){
 
    if(getShape() == nullptr){return false;}
 
    ivector1D dimLimit(3,2);
    //create internal shape using unique_ptr member.
    // unlink external shape eventually
    switch(getShapeType()){
        case ShapeType::CYLINDER :
            dimLimit[1] = 4;
            break;
        case ShapeType::SPHERE :
            dimLimit[1] = 5; dimLimit[2] = 3;
            break;
        default://CUBE //WEDGE
            break;
    }
 
    //check on dimensions and eventual closed loops on coordinates.
    m_nx = std::max(m_nx, dimLimit[0]-1);
    m_ny = std::max(m_ny, dimLimit[1]-1);
    m_nz = std::max(m_nz, dimLimit[2]-1);
 
    darray3E spanEff = getLocalSpan();
 
    reshapeNodalStructure();
 
    m_dx = spanEff[0]/m_nx;
    m_dy = spanEff[1]/m_ny;
    m_dz = spanEff[2]/m_nz;
    darray3E locOr = getShape()->getLocalOrigin();
    // get point distro;
    for (int i = 0; i < m_nx+1; i++) {m_xedge[i] = locOr[0] + ((double) i) * m_dx;}
    for (int i = 0; i < m_ny+1; i++) {m_yedge[i] = locOr[1] + ((double) i) * m_dy;}
    for (int i = 0; i < m_nz+1; i++) {m_zedge[i] = locOr[2] + ((double) i) * m_dz;}
    // get cell distro
    for (int i = 0; i < m_nx; i++) {m_xnode[i] = m_xedge[i] + 0.5 * m_dx;}
    for (int i = 0; i < m_ny; i++) {m_ynode[i] = m_yedge[i] + 0.5 * m_dy;}
    for (int i = 0; i < m_nz; i++) {m_znode[i] = m_zedge[i] + 0.5 * m_dz;}
 
    m_isBuild = true;
    return true;
};
 
 
void UStructMesh::execute(){
    build();
}
 
bool UStructMesh::isBuilt(){
    return(m_isBuild);
}
 
};