2D dynamic adaptive mesh refinement (AMR) using PABLO
2D dynamic adaptive mesh refinement (AMR) using PABLOThis example creates a 2D Octree mesh on the square domain [0,1]x[0,1].
The domain is adapted to track a set of moving bubbles with random initialization.
To run: ./PABLO_bubbles_2D
To see the result visit: PABLO website
/*---------------------------------------------------------------------------*\
 *
 *  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/>.
 *
\*---------------------------------------------------------------------------*/
 
#include "bitpit_PABLO.hpp"
 
#if BITPIT_ENABLE_MPI==1
#include "PABLO_userDataComm.hpp"
#include "PABLO_userDataLB.hpp"
#include <mpi.h>
#endif
 
using namespace std;
using namespace bitpit;
 
// =================================================================================== //
class bubble{
public:
    double c[2];
    double r;
};
 
void run()
{
    int iter = 0;
 
    PabloUniform pabloBB(2);
 
    pabloBB.setBalanceCodimension(1);
    uint32_t idx = 0;
    pabloBB.setBalance(idx,true);
 
    for (iter=1; iter<4; iter++){
        pabloBB.adaptGlobalRefine();
    }
 
#if BITPIT_ENABLE_MPI==1
    pabloBB.loadBalance();
#endif
 
    srand(0);
 
#if BITPIT_ENABLE_MPI==1
    int nb = 50;
#else
    int nb = 10;
#endif
    vector<bubble> BB;
    vector<bubble> BB0;
    vector<double> DY;
    vector<double> OM;
    vector<double> AA;
 
    for (int i=0; i<nb; i++){
        double randc[2];
        randc[0] = 0.8 * double(rand()) /  RAND_MAX + 0.1;
        randc[1] = double(rand()) /  RAND_MAX - 0.5;
        double randr = 0.1 * double(rand()) / RAND_MAX + 0.02;
        double dy = 0.005 + 0.05 * double(rand()) / RAND_MAX;
        double omega = 0.5 * double(rand()) / RAND_MAX;
        double aa = 0.15 * double(rand()) / RAND_MAX;
        bubble bb;
        bb.c[0] = randc[0];
        bb.c[1] = randc[1];
        bb.r = randr;
        BB.push_back(bb);
        BB0.push_back(bb);
        DY.push_back(dy);
        OM.push_back(omega);
        AA.push_back(aa);
    }
    double t0 = 0;
    double t = t0;
    double Dt = 0.5;
 
    int itstart = 1;
    int iterend = 200;
 
    for (iter=itstart; iter<iterend; iter++){
        if(pabloBB.getRank()==0) cout << "iter " << iter << endl;
        t += Dt;
 
        for (int i=0; i<nb; i++){
            BB[i].c[0] = BB0[i].c[0] + AA[i]*cos(OM[i]*t);
            BB[i].c[1] = BB[i].c[1]+ Dt*DY[i];
        }
 
        bool adapt = true;
        while (adapt){
            octantIterator it, itend = pabloBB.getInternalOctantsEnd();
            for (it=pabloBB.getInternalOctantsBegin(); it!=itend; ++it){
                bool inside = false;
                vector<array<double,3> > nodes = pabloBB.getNodes((*it));
                array<double,3> center = pabloBB.getCenter((*it));
                int ib = 0;
                while (!inside && ib<nb){
                    double xc = BB[ib].c[0];
                    double yc = BB[ib].c[1];
                    double radius = BB[ib].r;
                    for (int j=0; j<4; j++){
                        double x = nodes[j][0];
                        double y = nodes[j][1];
                        if ( ((!inside) &&
                                (pow((x-xc),2.0)+pow((y-yc),2.0) <= 1.25*pow(radius,2.0) &&
                                        pow((x-xc),2.0)+pow((y-yc),2.0) >= 0.75*pow(radius,2.0)))
                                        || ((!inside) && (pow((center[0]-xc),2.0)+pow((center[1]-yc),2.0) <= 1.25*pow(radius,2.0) &&
                                                pow((center[0]-xc),2.0)+pow((center[1]-yc),2.0) >= 0.75*pow(radius,2.0)))){
                            if (int(pabloBB.getLevel((*it))) < 9){
                                pabloBB.setMarker((*it),1);
                            }
                            else{
                                pabloBB.setMarker((*it),0);
                            }
                            inside = true;
                        }
                    }
                    ib++;
                }
                if (int(pabloBB.getLevel((*it))) > 0 && !inside){
                    pabloBB.setMarker((*it),5-pabloBB.getLevel((*it)));
                }
            }
 
            itend = pabloBB.getPboundOctantsEnd();
            for (it=pabloBB.getPboundOctantsBegin(); it!=itend; ++it){
                bool inside = false;
                vector<array<double,3> > nodes = pabloBB.getNodes((*it));
                array<double,3> center = pabloBB.getCenter((*it));
                int ib = 0;
                while (!inside && ib<nb){
                    double xc = BB[ib].c[0];
                    double yc = BB[ib].c[1];
                    double radius = BB[ib].r;
                    for (int j=0; j<4; j++){
                        double x = nodes[j][0];
                        double y = nodes[j][1];
                        if ( ((!inside) &&
                                (pow((x-xc),2.0)+pow((y-yc),2.0) <= 1.25*pow(radius,2.0) &&
                                        pow((x-xc),2.0)+pow((y-yc),2.0) >= 0.75*pow(radius,2.0)))
                                        || ((!inside) && (pow((center[0]-xc),2.0)+pow((center[1]-yc),2.0) <= 1.25*pow(radius,2.0) &&
                                                pow((center[0]-xc),2.0)+pow((center[1]-yc),2.0) >= 0.75*pow(radius,2.0)))){
                            if (pabloBB.getLevel((*it)) < 9){
                                pabloBB.setMarker((*it),1);
                            }
                            else{
                                pabloBB.setMarker((*it),0);
                            }
                            inside = true;
                        }
                    }
                    ib++;
                }
                if (pabloBB.getLevel((*it)) > 0 && !inside){
                    pabloBB.setMarker((*it),5-pabloBB.getLevel((*it)));
                }
            }
 
            adapt = pabloBB.adapt();
 
        }
 
#if BITPIT_ENABLE_MPI==1
        pabloBB.loadBalance();
#endif
 
        pabloBB.updateConnectivity();
        pabloBB.write("PabloBubble_iter"+to_string(static_cast<unsigned long long>(iter)));
 
    }
}
 
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
}