PABLO_example_00001.cpp

2D adaptive mesh refinement (AMR) using PABLO

2D adaptive mesh refinement (AMR) using PABLOThis example creates a 2D Octree mesh on the square domain [0,1]x[0,1].

The domain is refined globally one time, then refined iteratively using two different refinement criteria. At the end of the iterative refinement, one global refinement is performed again.

In the first criterion, each octant generated by a previous refinement iteration is marked for further refinement if its center is within a circle with a specified radius. Iterative refinement stops when a fixed number of iteration is reached.

In the second criterion, the refinement is performed using the same rules of the first criterion. However, this refinement is performed until the flag returned by the adapt method is true (for further details about adapt method, please visit documentation). Video below illustrates the iterative refinement. Each frame show the result of a refinement iteration starting from the ancestor quadrant up to the last global refinement.

The upper part of the domain is adapted by a refinement procedure using the first criterion, while the bottom part by a refinement procedure using the second criterion.

Moreover, in the right hand side of the domain, the 2:1 balancing is deactivated.

To run: ./PABLO_example_00001
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/>.
 *
\*---------------------------------------------------------------------------*/
 
#if BITPIT_ENABLE_MPI==1
#include <mpi.h>
#endif
 
#include "bitpit_PABLO.hpp"
 
using namespace std;
using namespace bitpit;
 
// =================================================================================== //
// =================================================================================== //
 
void run()
{
    PabloUniform pablo1(2);
 
    pablo1.computeConnectivity();
    pablo1.write("pablo00001_iter0");
 
    pablo1.adaptGlobalRefine();
    pablo1.updateConnectivity();
    pablo1.write("pablo00001_iter1");
 
    double xc, yc;
    xc = yc = 0.5;
 
    pablo1.setBalanceCodimension(1);
 
    uint32_t nocts = pablo1.getNumOctants();
    for (unsigned int i=0; i<nocts; i++){
        array<double,3> center = pablo1.getCenter(i);
        double x = center[0];
        if (x>xc)
            pablo1.setBalance(i,false);
    }
 
    double radius = 0.4;
 
    int nref1 = 6;
    for (int iter=0; iter<nref1; iter++){
        nocts = pablo1.getNumOctants();
        for (unsigned int i=0; i<nocts; i++){
            Octant *oct = pablo1.getOctant(i);
            array<double,3> center = pablo1.getCenter(oct);
            double x = center[0];
            double y = center[1];
 
            if ((pow((x-xc),2.0)+pow((y-yc),2.0) < pow(radius,2.0)) &&
                    (y<yc)){
                pablo1.setMarker(oct, 1);
            }
        }
        pablo1.adapt();
        pablo1.updateConnectivity();
        pablo1.write("pablo00001_iter"+to_string(static_cast<unsigned long long>(iter+2)));
    }
 
    int nref2 = 5;
    int iter = 0;
    bool done = true;
    while(iter<=nref2){
        done = true;
        while(done)
        {
            nocts = pablo1.getNumOctants();
            for (unsigned int i=0; i<nocts; i++){
                array<double,3> center = pablo1.getCenter(i);
                double x = center[0];
                double y = center[1];
                if ((pow((x-xc),2.0)+pow((y-yc),2.0) < pow(radius,2.0)) &&
                        (y>yc) && iter<=nref2 && pablo1.getLevel(i)<=iter+1){
 
                    pablo1.setMarker(i, 1);
                }
            }
            done = pablo1.adapt();
            pablo1.updateConnectivity();
            pablo1.write("pablo00001_iter"+to_string(static_cast<unsigned long long>(iter+nref1+2)));
        }
        iter++;
    }
    pablo1.adaptGlobalRefine();
    pablo1.updateConnectivity();
    pablo1.write("pablo00001_iter"+to_string(static_cast<unsigned long long>(iter+nref1+3)));
}
 
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
}