PABLO_example_00006.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/>.
 *
\*---------------------------------------------------------------------------*/
 
#if BITPIT_ENABLE_MPI==1
#include <mpi.h>
#endif
 
#include "bitpit_PABLO.hpp"
 
#if BITPIT_ENABLE_MPI==1
#include "PABLO_userDataComm.hpp"
#endif
 
using namespace std;
using namespace bitpit;
 
// =================================================================================== //
// =================================================================================== //
 
class Data{
public:
    vector<double> doubleData;
    vector<float> floatData;
    Data(uint32_t nocts): doubleData(nocts,0.0), floatData(nocts,0.0){};
    Data(Data& rhs){
        this->doubleData = rhs.doubleData;
        this->floatData = rhs.floatData;
    }
};
 
void run()
{
    int iter = 0;
    int dim = 2;
 
    PabloUniform pablo6(2);
 
    for (iter=1; iter<5; iter++){
        pablo6.adaptGlobalRefine();
    }
 
#if BITPIT_ENABLE_MPI==1
    pablo6.loadBalance();
#endif
 
    double xc, yc;
    xc = yc = 0.5;
    double radius = 0.25;
 
    uint32_t nocts = pablo6.getNumOctants();
    uint32_t nghosts = pablo6.getNumGhosts();
    //vector<double> oct_data(nocts, 0.0), ghost_data(nghosts, 0.0);
    Data octdata(nocts), ghostdata(nghosts);
 
 
    for (unsigned int i=0; i<nocts; i++){
        vector<array<double,3> > nodes = pablo6.getNodes(i);
        for (int j=0; j<4; j++){
            double x = nodes[j][0];
            double y = nodes[j][1];
            if ((pow((x-xc),2.0)+pow((y-yc),2.0) <= pow(radius,2.0))){
                //oct_data[i] = 1.0;
                octdata.doubleData[i] = 1.0;
                octdata.floatData[i] = 1.0;
            }
        }
    }
 
    for (unsigned int i=0; i<nghosts; i++){
        Octant *oct = pablo6.getGhostOctant(i);
        vector<array<double,3> > nodes = pablo6.getNodes(oct);
        for (int j=0; j<4; j++){
            double x = nodes[j][0];
            double y = nodes[j][1];
            if ((pow((x-xc),2.0)+pow((y-yc),2.0) <= pow(radius,2.0))){
                //ghost_data[i] = 1.0;
                ghostdata.doubleData[i] = 1.0;
                ghostdata.floatData[i] = 1.0;
            }
        }
    }
 
    iter = 0;
    pablo6.updateConnectivity();
    pablo6.writeTest("pablo00006_double_iter"+to_string(static_cast<unsigned long long>(iter)), octdata.doubleData);
 
    int start = iter + 1;
    for (iter=start; iter<start+25; iter++){
        //vector<double> oct_data_smooth(nocts, 0.0);
        Data octdatasmooth(nocts);
        vector<uint32_t> neigh, neigh_t;
        vector<bool> isghost, isghost_t;
        uint8_t iface, nfaces, codim;
        for (unsigned int i=0; i<nocts; i++){
            neigh.clear();
            isghost.clear();
 
            for (codim=1; codim<dim+1; codim++){
                if (codim == 1){
                    nfaces = 4;
                }
                else if (codim == 2){
                    nfaces = 4;
                }
                for (iface=0; iface<nfaces; iface++){
                    pablo6.findNeighbours(i,iface,codim,neigh_t,isghost_t);
                    neigh.insert(neigh.end(), neigh_t.begin(), neigh_t.end());
                    isghost.insert(isghost.end(), isghost_t.begin(), isghost_t.end());
                }
            }
 
            //oct_data_smooth[i] = oct_data[i]/(neigh.size()+1);
            octdatasmooth.doubleData[i] = octdata.doubleData[i]/(neigh.size()+1);
            octdatasmooth.floatData[i] = octdata.floatData[i]/(neigh.size()+1);
            for (unsigned int j=0; j<neigh.size(); j++){
                if (isghost[j]){
                    //oct_data_smooth[i] += ghost_data[neigh[j]]/(neigh.size()+1);
                    octdatasmooth.doubleData[i] += ghostdata.doubleData[neigh[j]]/(neigh.size()+1);
                    octdatasmooth.floatData[i] += ghostdata.floatData[neigh[j]]/(neigh.size()+1);
                }
                else{
                    //oct_data_smooth[i] += oct_data[neigh[j]]/(neigh.size()+1);
                    octdatasmooth.doubleData[i] += octdata.doubleData[neigh[j]]/(neigh.size()+1);
                    octdatasmooth.floatData[i] += octdata.floatData[neigh[j]]/(neigh.size()+1);
                }
            }
        }
 
        pablo6.updateConnectivity();
        pablo6.writeTest("pablo00006_iter"+to_string(static_cast<unsigned long long>(iter)), octdatasmooth.doubleData);
 
#if BITPIT_ENABLE_MPI==1
        UserDataComm<Data> data_comm(octdatasmooth, ghostdata);
        pablo6.communicate(data_comm);
 
#endif
        octdata.doubleData = octdatasmooth.doubleData;
        octdata.floatData = octdatasmooth.floatData;
    }
}
 
int main(int argc, char *argv[])
{
#if BITPIT_ENABLE_MPI==1
    MPI_Init(&argc,&argv);
#else
    BITPIT_UNUSED(argc);
    BITPIT_UNUSED(argv);
#endif
 
    int nProcs;
    int rank;
#if BITPIT_ENABLE_MPI==1
    MPI_Comm_size(MPI_COMM_WORLD, &nProcs);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
#else
    nProcs = 1;
    rank   = 0;
#endif
 
    // Initialize the logger
    log::manager().initialize(log::MODE_SEPARATE, false, nProcs, rank);
    log::cout() << log::fileVerbosity(log::LEVEL_INFO);
    log::cout() << log::disableConsole();
 
    // Run the example
    try {
        run();
    } catch (const std::exception &exception) {
        log::cout() << exception.what();
        exit(1);
    }
 
#if BITPIT_ENABLE_MPI==1
    MPI_Finalize();
#endif
}