PabloUniform.cpp

/*---------------------------------------------------------------------------*\
 *
 *  bitpit
 *
 *  Copyright (C) 2015-2021 OPTIMAD engineering Srl
 *
 *  -------------------------------------------------------------------------
 *  License
 *  This file is part of bitpit.
 *
 *  bitpit 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.
 *
 *  bitpit 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 bitpit. If not, see <http://www.gnu.org/licenses/>.
 *
 \*---------------------------------------------------------------------------*/
 
// =================================================================================== //
// INCLUDES                                                                            //
// =================================================================================== //
#include "bitpit_operators.hpp"
 
#include "PabloUniform.hpp"
 
namespace bitpit {
 
    // =================================================================================== //
    // NAME SPACES                                                                         //
    // =================================================================================== //
    using namespace std;
 
    // =================================================================================== //
    // CLASS IMPLEMENTATION                                                                    //
    // =================================================================================== //
 
    // =================================================================================== //
    // CONSTRUCTORS AND OPERATORS
    // =================================================================================== //
#if BITPIT_ENABLE_MPI==1
    PabloUniform::PabloUniform(const std::string &logfile, MPI_Comm comm):ParaTree(logfile,comm){
#else
    PabloUniform::PabloUniform(const std::string &logfile):ParaTree(logfile){
#endif
        __reset();
    }
 
#if BITPIT_ENABLE_MPI==1
    PabloUniform::PabloUniform(uint8_t dim, const std::string &logfile, MPI_Comm comm):ParaTree(dim,logfile,comm){
#else
    PabloUniform::PabloUniform(uint8_t dim, const std::string &logfile):ParaTree(dim,logfile){
#endif
        __reset();
    };
 
#if BITPIT_ENABLE_MPI==1
    PabloUniform::PabloUniform(double X, double Y, double Z, double L, uint8_t dim, const std::string &logfile, MPI_Comm comm):ParaTree(dim,logfile,comm){
#else
    PabloUniform::PabloUniform(double X, double Y, double Z, double L, uint8_t dim, const std::string &logfile):ParaTree(dim,logfile){
#endif
        __reset();
 
        setOrigin({{X, Y, Z}});
        setL(L);
    };
 
    // =================================================================================== //
    // METHODS
    // =================================================================================== //
 
    void
    PabloUniform::reset(){
        ParaTree::reset();
        __reset();
    }
 
    void
    PabloUniform::__reset(){
        setOrigin({{0,0,0}});
        setL(1.);
    }
 
    int
    PabloUniform::getDumpVersion() const
    {
        const int DUMP_VERSION = 1;
 
        return (DUMP_VERSION + ParaTree::getDumpVersion());
    }
 
    void
    PabloUniform::dump(std::ostream &stream, bool full)
    {
        ParaTree::dump(stream, full);
 
        utils::binary::write(stream, m_origin[0]);
        utils::binary::write(stream, m_origin[1]);
        utils::binary::write(stream, m_origin[2]);
        utils::binary::write(stream, m_L);
    }
 
    void
    PabloUniform::restore(std::istream &stream)
    {
        ParaTree::restore(stream);
 
        std::array<double, 3> origin = {{0., 0., 0.}};
        utils::binary::read(stream, origin[0]);
        utils::binary::read(stream, origin[1]);
        utils::binary::read(stream, origin[2]);
        setOrigin(origin);
 
        double L;
        utils::binary::read(stream, L);
        setL(L);
    }
 
    // =================================================================================== //
    // BASIC GET/SET METHODS                                                               //
    // =================================================================================== //
    darray3
    PabloUniform::getOrigin() const {
        return m_origin;
    };
 
    double
    PabloUniform::getX0() const {
        return m_origin[0];
    };
 
    double
    PabloUniform::getY0() const {
        return m_origin[1];
    };
 
    double
    PabloUniform::getZ0() const {
        return m_origin[2];
    };
 
    double
    PabloUniform::getL() const {
        return m_L;
    };
 
    void
    PabloUniform::setL(double L){
        m_L      = L;
        m_area   = uipow(L, getDim() - 1);
        m_volume = uipow(L, getDim());
    };
 
    void
    PabloUniform::setOrigin(const darray3 &origin){
        m_origin = origin;
    };
 
    double
    PabloUniform::levelToSize(uint8_t level) const {
        double size = ParaTree::levelToSize(level);
        return m_L *size;
    }
 
    // =================================================================================== //
    // INDEX BASED METHODS                                                                 //
    // =================================================================================== //
    darray3
    PabloUniform::getCoordinates(uint32_t idx) const {
        darray3 coords  = {{0., 0., 0.}};
        darray3 coords_ = ParaTree::getCoordinates(idx);
        for (int i=0; i<3; i++){
            coords[i] = m_origin[i] + m_L * coords_[i];
        }
        return coords;
    };
 
    double
    PabloUniform::getX(uint32_t idx) const {
        double X, X_;
        X_ = ParaTree::getX(idx);
        X = m_origin[0] + m_L * X_;
        return X;
    };
 
    double
    PabloUniform::getY(uint32_t idx) const {
        double Y, Y_;
        Y_ = ParaTree::getY(idx);
        Y = m_origin[1] + m_L * Y_;
        return Y;
    };
 
    double
    PabloUniform::getZ(uint32_t idx) const {
        double Z, Z_;
        Z_ = ParaTree::getZ(idx);
        Z = m_origin[2] + m_L * Z_;
        return Z;
    };
 
    double
    PabloUniform::getSize(uint32_t idx) const {
        return m_L * ParaTree::getSize(idx);
    };
 
    double
    PabloUniform::getArea(uint32_t idx) const {
        return m_area * ParaTree::getArea(idx);
    };
 
    double
    PabloUniform::getVolume(uint32_t idx) const {
        return m_volume * ParaTree::getVolume(idx);
    };
 
    void
    PabloUniform::getCenter(uint32_t idx, darray3& center) const {
        darray3 center_ = ParaTree::getCenter(idx);
        for (int i=0; i<3; i++){
            center[i] = m_origin[i] + m_L * center_[i];
        }
    };
 
    darray3
    PabloUniform::getCenter(uint32_t idx) const {
        darray3 center  = {{0., 0., 0.}};
        darray3 center_ = ParaTree::getCenter(idx);
        for (int i=0; i<3; i++){
            center[i] = m_origin[i] + m_L * center_[i];
        }
        return center;
    };
 
    void
    PabloUniform::getFaceCenter(uint32_t idx, uint8_t iface, darray3& center) const {
        darray3 center_ = ParaTree::getFaceCenter(idx, iface);
        for (int i=0; i<3; i++){
            center[i] = m_origin[i] + m_L * center_[i];
        }
    };
 
    darray3
    PabloUniform::getFaceCenter(uint32_t idx, uint8_t iface) const {
        darray3 center  = {{0., 0., 0.}};
        darray3 center_ = ParaTree::getFaceCenter(idx, iface);
        for (int i=0; i<3; i++){
            center[i] = m_origin[i] + m_L * center_[i];
        }
        return center;
    };
 
    darray3
    PabloUniform::getNode(uint32_t idx, uint8_t inode) const {
        darray3 node  = {{0., 0., 0.}};
        darray3 node_ = ParaTree::getNode(idx, inode);
        for (int i=0; i<3; i++){
            node[i] = m_origin[i] + m_L * node_[i];
        }
        return node;
    };
 
    void
    PabloUniform::getNode(uint32_t idx, uint8_t inode, darray3& node) const {
        darray3 node_ = ParaTree::getNode(idx, inode);
        for (int i=0; i<3; i++){
            node[i] = m_origin[i] + m_L * node_[i];
        }
    };
 
    void
    PabloUniform::getNodes(uint32_t idx, darr3vector & nodes) const {
        darray3vector nodes_ = ParaTree::getNodes(idx);
        nodes.resize(ParaTree::getNnodes());
        for (int j=0; j<ParaTree::getNnodes(); j++){
            for (int i=0; i<3; i++){
                nodes[j][i] = m_origin[i] + m_L * nodes_[j][i];
            }
        }
    };
 
    darr3vector
    PabloUniform::getNodes(uint32_t idx) const {
        darray3vector nodes, nodes_ = ParaTree::getNodes(idx);
        nodes.resize(ParaTree::getNnodes());
        for (int j=0; j<ParaTree::getNnodes(); j++){
            for (int i=0; i<3; i++){
                nodes[j][i] = m_origin[i] + m_L * nodes_[j][i];
            }
        }
        return nodes;
    };
 
    void
    PabloUniform::getNormal(uint32_t idx, uint8_t iface, darray3 & normal) const {
        ParaTree::getNormal(idx, iface, normal);
    }
 
    darray3
    PabloUniform::getNormal(uint32_t idx, uint8_t iface) const {
        return ParaTree::getNormal(idx, iface);
    }
 
    // =================================================================================== //
    // POINTER BASED METHODS                                                               //
    // =================================================================================== //
    darray3
    PabloUniform::getCoordinates(const Octant* oct) const {
        darray3 coords  = {{0., 0., 0.}};
        darray3 coords_ = ParaTree::getCoordinates(oct);
        for (int i=0; i<3; i++){
            coords[i] = m_origin[i] + m_L * coords_[i];
        }
        return coords;
    };
 
    double
    PabloUniform::getX(const Octant* oct) const {
        double X, X_;
        X_ = ParaTree::getX(oct);
        X = m_origin[0] + m_L * X_;
        return X;
    };
 
    double
    PabloUniform::getY(const Octant* oct) const {
        double Y, Y_;
        Y_ = ParaTree::getY(oct);
        Y = m_origin[1] + m_L * Y_;
        return Y;
    };
 
    double
    PabloUniform::getZ(const Octant* oct) const {
        double Z, Z_;
        Z_ = ParaTree::getZ(oct);
        Z = m_origin[2] + m_L * Z_;
        return Z;
    };
 
    double
    PabloUniform::getSize(const Octant* oct) const {
        return m_L * ParaTree::getSize(oct);
    };
 
    double
    PabloUniform::getArea(const Octant* oct) const {
        return m_area * ParaTree::getArea(oct);
    };
 
    double
    PabloUniform::getVolume(const Octant* oct) const {
        return m_volume * ParaTree::getVolume(oct);
    };
 
    void
    PabloUniform::getCenter(const Octant* oct, darray3& center) const {
        darray3 center_ = ParaTree::getCenter(oct);
        for (int i=0; i<3; i++){
            center[i] = m_origin[i] + m_L * center_[i];
        }
    };
 
    darray3
    PabloUniform::getCenter(const Octant* oct) const {
        darray3 center  = {{0., 0., 0.}};
        darray3 center_ = ParaTree::getCenter(oct);
        for (int i=0; i<3; i++){
            center[i] = m_origin[i] + m_L * center_[i];
        }
        return center;
    };
 
    void
    PabloUniform::getFaceCenter(const Octant* oct, uint8_t iface, darray3& center) const {
        darray3 center_ = ParaTree::getFaceCenter(oct, iface);
        for (int i=0; i<3; i++){
            center[i] = m_origin[i] + m_L * center_[i];
        }
    };
 
    darray3
    PabloUniform::getFaceCenter(const Octant* oct, uint8_t iface) const {
        darray3 center  = {{0., 0., 0.}};
        darray3 center_ = ParaTree::getFaceCenter(oct, iface);
        for (int i=0; i<3; i++){
            center[i] = m_origin[i] + m_L * center_[i];
        }
        return center;
    };
 
    darray3
    PabloUniform::getNode(const Octant* oct, uint8_t inode) const {
        darray3 node  = {{0., 0., 0.}};
        darray3 node_ = ParaTree::getNode(oct, inode);
        for (int i=0; i<3; i++){
            node[i] = m_origin[i] + m_L * node_[i];
        }
        return node;
    };
 
    void
    PabloUniform::getNode(const Octant* oct, uint8_t inode, darray3& node) const {
        darray3 node_ = ParaTree::getNode(oct, inode);
        for (int i=0; i<3; i++){
            node[i] = m_origin[i] + m_L * node_[i];
        }
    };
 
    void
    PabloUniform::getNodes(const Octant* oct, darr3vector & nodes) const {
        darray3vector nodes_ = ParaTree::getNodes(oct);
        nodes.resize(ParaTree::getNnodes());
        for (int j=0; j<ParaTree::getNnodes(); j++){
            for (int i=0; i<3; i++){
                nodes[j][i] = m_origin[i] + m_L * nodes_[j][i];
            }
        }
    };
 
    darr3vector
    PabloUniform::getNodes(const Octant* oct) const {
        darray3vector nodes  = {{{{0., 0., 0.}}, {{0., 0., 0.}}, {{0., 0., 0.}}}};
        darray3vector nodes_ = ParaTree::getNodes(oct);
        nodes.resize(ParaTree::getNnodes());
        for (int j=0; j<ParaTree::getNnodes(); j++){
            for (int i=0; i<3; i++){
                nodes[j][i] = m_origin[i] + m_L * nodes_[j][i];
            }
        }
        return nodes;
    };
 
    void
    PabloUniform::getNormal(const Octant* oct, uint8_t iface, darray3 & normal) const {
        ParaTree::getNormal(oct, iface, normal);
    }
 
    darray3
    PabloUniform::getNormal(const Octant* oct, uint8_t iface) const {
        return ParaTree::getNormal(oct, iface);
    }
 
    // =================================================================================== //
    // LOCAL TREE GET/SET METHODS                                                          //
    // =================================================================================== //
    double
    PabloUniform::getLocalMaxSize() const {
        return m_L * ParaTree::getLocalMaxSize();
    };
 
    double
    PabloUniform::getLocalMinSize() const {
        return m_L * ParaTree::getLocalMinSize();
    };
 
 
    void
    PabloUniform::getBoundingBox(darray3 & P0, darray3 & P1) const {
        // If there are no octants the bounding box is empty
        uint32_t nocts = ParaTree::getNumOctants();
        if (nocts == 0) {
            P0 = getOrigin();
            P1 = P0;
 
            return;
        }
 
        // If the octree is serial we can evaluate the bounding box easily
        // otherwise we need to scan all the octants
        if (getSerial()) {
            P0 = getOrigin();
            P1 = P0;
            for (int i=0; i<ParaTree::getDim(); i++){
                P1[i] += getL();
            }
 
            return;
        }
 
        // If the octree is parallel we need to scan all the octants
        darray3     cnode0, cnode1;
 
        uint32_t    id = 0;
        uint8_t     nnodes = ParaTree::getNnodes();
 
        P0 = getNode(id, 0);
        P1 = getNode(nocts-1, nnodes-1);
 
        for (id=0; id<nocts; id++){
            cnode0 = getNode(id, 0);
            cnode1 = getNode(id, nnodes-1);
            for (int i=0; i<ParaTree::getDim(); i++){
                P0[i] = min(P0[i], cnode0[i]);
                P1[i] = max(P1[i], cnode1[i]);
            }
        }
    };
 
 
    // =================================================================================== //
    // INTERSECTION GET/SET METHODS                                                        //
    // =================================================================================== //
    double
    PabloUniform::getSize(const Intersection* inter) const {
        return m_L * ParaTree::getSize(inter);
    };
 
    double
    PabloUniform::getArea(const Intersection* inter) const {
        return m_area * ParaTree::getArea(inter);
    };
 
    darray3
    PabloUniform::getCenter(const Intersection* inter) const {
        darray3 center = ParaTree::getCenter(inter);
        for (int i=0; i<3; i++){
            center[i] = m_origin[i] + m_L * center[i];
        }
        return center;
    }
 
    darr3vector
    PabloUniform::getNodes(const Intersection* inter) const {
        darr3vector nodes, nodes_ = ParaTree::getNodes(inter);
        nodes.resize(ParaTree::getNnodesperface());
        for (int j=0; j<ParaTree::getNnodesperface(); j++){
            for (int i=0; i<3; i++){
                nodes[j][i] = m_origin[i] + m_L * nodes_[j][i];
            }
        }
        return nodes;
    }
 
    darray3
    PabloUniform::getNormal(const Intersection* inter) const {
        return ParaTree::getNormal(inter);
    }
 
    // =================================================================================== //
    // OTHER OCTANT BASED METHODS                                                          //
    // =================================================================================== //
    Octant* PabloUniform::getPointOwner(darray3 point){
        for (int i=0; i<3; i++){
            point[i] = (point[i] - m_origin[i])/m_L;
        }
        return ParaTree::getPointOwner(point);
    };
 
    Octant* PabloUniform::getPointOwner(darray3 point, bool & isghost){
        for (int i=0; i<3; i++){
            point[i] = (point[i] - m_origin[i])/m_L;
        }
        return ParaTree::getPointOwner(point,isghost);
    };
 
 
    uint32_t
    PabloUniform::getPointOwnerIdx(darray3 point) const {
        for (int i=0; i<3; i++){
            point[i] = (point[i] - m_origin[i])/m_L;
        }
        return ParaTree::getPointOwnerIdx(point);
    };
    
    uint32_t
    PabloUniform::getPointOwnerIdx(darray3 point, bool & isghost) const {
        for (int i=0; i<3; i++){
            point[i] = (point[i] - m_origin[i])/m_L;
        }
        return ParaTree::getPointOwnerIdx(point,isghost);
    };
    
    int PabloUniform::getPointOwnerRank(darray3 point){
        for (int i=0; i<3; i++){
            point[i] = (point[i] - m_origin[i])/m_L;
        }
        return ParaTree::getPointOwnerRank(point);
    };
    
    // =================================================================================== //
    // OTHER PARATREE BASED METHODS                                                        //
    // =================================================================================== //
    darray3
    PabloUniform::getNodeCoordinates(uint32_t inode) const {
        darray3 node = ParaTree::getNodeCoordinates(inode);
        for (int i=0; i<3; i++){
            node[i] = m_origin[i] + m_L * node[i];
        }
        return node;
    }
 
    // =================================================================================== //
    // TESTING OUTPUT METHODS                                                              //
    // =================================================================================== //
    void
    PabloUniform::write(const std::string &filename) {
 
        if (getConnectivity().size() == 0) {
            computeConnectivity();
        }
 
        stringstream name;
        name << "s" << std::setfill('0') << std::setw(4) << getNproc() << "-p" << std::setfill('0') << std::setw(4) << getRank() << "-" << filename << ".vtu";
 
        ofstream out(name.str().c_str());
        if(!out.is_open()){
            stringstream ss;
            ss << filename << "*.vtu cannot be opened and it won't be written." << endl;
            getLog() << ss.str();
            return;
        }
        int nofNodes = getNumNodes();
        int nofOctants = getNumOctants();
        int nofGhosts = getNumGhosts();
        int nofAll = nofGhosts + nofOctants;
        out << "<?xml version=\"1.0\"?>" << endl
            << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">" << endl
            << "  <UnstructuredGrid>" << endl
            << "    <Piece NumberOfCells=\"" << getConnectivity().size() + getGhostConnectivity().size() << "\" NumberOfPoints=\"" << getNumNodes() << "\">" << endl;
        out << "      <Points>" << endl
            << "        <DataArray type=\"Float64\" Name=\"Coordinates\" NumberOfComponents=\""<< 3 <<"\" format=\"ascii\">" << endl
            << "          " << std::fixed;
        for(int i = 0; i < nofNodes; i++)
            {
                const std::array<double,3> & nodeCoordinates = getNodeCoordinates(i);
                for(int j = 0; j < 3; ++j){
                    out << std::setprecision(6) << nodeCoordinates[j] << " ";
                }
                if((i+1)%4==0 && i!=nofNodes-1)
                    out << endl << "          ";
            }
        out << endl << "        </DataArray>" << endl
            << "      </Points>" << endl
            << "      <Cells>" << endl
            << "        <DataArray type=\"UInt64\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"ascii\">" << endl
            << "          ";
        for(int i = 0; i < nofOctants; i++)
            {
                for(int j = 0; j < getNnodes(); j++)
                    {
                        int jj = j;
                        if (getDim()==2){
                            if (j<2){
                                jj = j;
                            }
                            else if(j==2){
                                jj = 3;
                            }
                            else if(j==3){
                                jj = 2;
                            }
                        }
                        out << getConnectivity()[i][jj] << " ";
                    }
                if((i+1)%3==0 && i!=nofOctants-1)
                    out << endl << "          ";
            }
        for(int i = 0; i < nofGhosts; i++)
            {
                for(int j = 0; j < getNnodes(); j++)
                    {
                        int jj = j;
                        if (getDim()==2){
                            if (j<2){
                                jj = j;
                            }
                            else if(j==2){
                                jj = 3;
                            }
                            else if(j==3){
                                jj = 2;
                            }
                        }
                        out << getGhostConnectivity()[i][jj] << " ";
                    }
                if((i+1)%3==0 && i!=nofGhosts-1)
                    out << endl << "          ";
            }
        out << endl << "        </DataArray>" << endl
            << "        <DataArray type=\"UInt64\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"ascii\">" << endl
            << "          ";
        for(int i = 0; i < nofAll; i++)
            {
                out << (i+1)*getNnodes() << " ";
                if((i+1)%12==0 && i!=nofAll-1)
                    out << endl << "          ";
            }
        out << endl << "        </DataArray>" << endl
            << "        <DataArray type=\"UInt8\" Name=\"types\" NumberOfComponents=\"1\" format=\"ascii\">" << endl
            << "          ";
        for(int i = 0; i < nofAll; i++)
            {
                int type;
                type = 5 + (getDim()*2);
                out << type << " ";
                if((i+1)%12==0 && i!=nofAll-1)
                    out << endl << "          ";
            }
        out << endl << "        </DataArray>" << endl
            << "      </Cells>" << endl
            << "    </Piece>" << endl
            << "  </UnstructuredGrid>" << endl
            << "</VTKFile>" << endl;
 
 
        if(getRank() == 0){
            name.str("");
            name << "s" << std::setfill('0') << std::setw(4) << getNproc() << "-" << filename << ".pvtu";
            ofstream pout(name.str().c_str());
            if(!pout.is_open()){
                stringstream ss;
                ss << filename << "*.pvtu cannot be opened and it won't be written." << endl;
                getLog() << ss.str();
                return;
            }
 
            pout << "<?xml version=\"1.0\"?>" << endl
                 << "<VTKFile type=\"PUnstructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">" << endl
                 << "  <PUnstructuredGrid GhostLevel=\"0\">" << endl
                 << "    <PPointData>" << endl
                 << "    </PPointData>" << endl
                 << "    <PCellData Scalars=\"\">" << endl;
            pout << "    </PCellData>" << endl
                 << "    <PPoints>" << endl
                 << "      <PDataArray type=\"Float64\" Name=\"Coordinates\" NumberOfComponents=\"3\"/>" << endl
                 << "    </PPoints>" << endl;
            for(int i = 0; i < getNproc(); i++)
                pout << "    <Piece Source=\"s" << std::setw(4) << std::setfill('0') << getNproc() << "-p" << std::setw(4) << std::setfill('0') << i << "-" << filename << ".vtu\"/>" << endl;
            pout << "  </PUnstructuredGrid>" << endl
                 << "</VTKFile>";
 
            pout.close();
 
        }
#if BITPIT_ENABLE_MPI==1
        if (isCommSet()) {
            MPI_Barrier(getComm());
        }
#endif
 
    }
 
    void
    PabloUniform::writeTest(const std::string &filename, vector<double> data) {
 
        if (getConnectivity().size() == 0) {
            computeConnectivity();
        }
 
        stringstream name;
        name << "s" << std::setfill('0') << std::setw(4) << getNproc() << "-p" << std::setfill('0') << std::setw(4) << getRank() << "-" << filename << ".vtu";
 
        ofstream out(name.str().c_str());
        if(!out.is_open()){
            stringstream ss;
            ss << filename << "*.vtu cannot be opened and it won't be written.";
            getLog() << ss.str();
            return;
        }
        int nofNodes = getNumNodes();
        int nofOctants = getNumOctants();
        int nofAll = nofOctants;
        out << "<?xml version=\"1.0\"?>" << endl
            << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">" << endl
            << "  <UnstructuredGrid>" << endl
            << "    <Piece NumberOfCells=\"" << getNumOctants() << "\" NumberOfPoints=\"" << getNumNodes() << "\">" << endl;
        out << "      <CellData Scalars=\"Data\">" << endl;
        out << "      <DataArray type=\"Float64\" Name=\"Data\" NumberOfComponents=\"1\" format=\"ascii\">" << endl
            << "          " << std::fixed;
        int ndata = getNumOctants();
        for(int i = 0; i < ndata; i++)
            {
                out << std::setprecision(6) << data[i] << " ";
                if((i+1)%4==0 && i!=ndata-1)
                    out << endl << "          ";
            }
        out << endl << "        </DataArray>" << endl
            << "      </CellData>" << endl
            << "      <Points>" << endl
            << "        <DataArray type=\"Float64\" Name=\"Coordinates\" NumberOfComponents=\""<< 3 <<"\" format=\"ascii\">" << endl
            << "          " << std::fixed;
        for(int i = 0; i < nofNodes; i++)
            {
            const std::array<double,3> & nodeCoordinates = getNodeCoordinates(i);
            for(int j = 0; j < 3; ++j){
                out << std::setprecision(6) << nodeCoordinates[j] << " ";
            }
                if((i+1)%4==0 && i!=nofNodes-1)
                    out << endl << "          ";
            }
        out << endl << "        </DataArray>" << endl
            << "      </Points>" << endl
            << "      <Cells>" << endl
            << "        <DataArray type=\"UInt64\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"ascii\">" << endl
            << "          ";
        for(int i = 0; i < nofOctants; i++)
            {
                for(int j = 0; j < getNnodes(); j++)
                    {
                        int jj = j;
                        if (getDim()==2){
                            if (j<2){
                                jj = j;
                            }
                            else if(j==2){
                                jj = 3;
                            }
                            else if(j==3){
                                jj = 2;
                            }
                        }
                        out << getConnectivity()[i][jj] << " ";
                    }
                if((i+1)%3==0 && i!=nofOctants-1)
                    out << endl << "          ";
            }
        out << endl << "        </DataArray>" << endl
            << "        <DataArray type=\"UInt64\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"ascii\">" << endl
            << "          ";
        for(int i = 0; i < nofAll; i++)
            {
                out << (i+1)*getNnodes() << " ";
                if((i+1)%12==0 && i!=nofAll-1)
                    out << endl << "          ";
            }
        out << endl << "        </DataArray>" << endl
            << "        <DataArray type=\"UInt8\" Name=\"types\" NumberOfComponents=\"1\" format=\"ascii\">" << endl
            << "          ";
        for(int i = 0; i < nofAll; i++)
            {
                int type;
                type = 5 + (getDim()*2);
                out << type << " ";
                if((i+1)%12==0 && i!=nofAll-1)
                    out << endl << "          ";
            }
        out << endl << "        </DataArray>" << endl
            << "      </Cells>" << endl
            << "    </Piece>" << endl
            << "  </UnstructuredGrid>" << endl
            << "</VTKFile>" << endl;
 
 
        if(getRank() == 0){
            name.str("");
            name << "s" << std::setfill('0') << std::setw(4) << getNproc() << "-" << filename << ".pvtu";
            ofstream pout(name.str().c_str());
            if(!pout.is_open()){
                stringstream ss;
                ss << filename << "*.pvtu cannot be opened and it won't be written." << endl;
                getLog() << ss.str();
                return;
            }
 
            pout << "<?xml version=\"1.0\"?>" << endl
                 << "<VTKFile type=\"PUnstructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">" << endl
                 << "  <PUnstructuredGrid GhostLevel=\"0\">" << endl
                 << "    <PPointData>" << endl
                 << "    </PPointData>" << endl
                 << "    <PCellData Scalars=\"Data\">" << endl
                 << "      <PDataArray type=\"Float64\" Name=\"Data\" NumberOfComponents=\"1\"/>" << endl
                 << "    </PCellData>" << endl
                 << "    <PPoints>" << endl
                 << "      <PDataArray type=\"Float64\" Name=\"Coordinates\" NumberOfComponents=\"3\"/>" << endl
                 << "    </PPoints>" << endl;
            for(int i = 0; i < getNproc(); i++)
                pout << "    <Piece Source=\"s" << std::setw(4) << std::setfill('0') << getNproc() << "-p" << std::setw(4) << std::setfill('0') << i << "-" << filename << ".vtu\"/>" << endl;
            pout << "  </PUnstructuredGrid>" << endl
                 << "</VTKFile>";
 
            pout.close();
 
        }
#if BITPIT_ENABLE_MPI==1
        MPI_Barrier(getComm());
#endif
 
    }
 
}