PABLO_example_00004.cpp

2D adaptive mesh refinement (AMR) with data using PABLO

2D adaptive mesh refinement (AMR) with data using PABLOThis example shows how to use PABLO for adaptive mesh refinement with data associated to the mesh.

A set of data (value 0 assigned to all quadrants outside a circle, and value equal to the distance from circle's center assigned to all quadrants within the circle) is linked to a uniform quadtree of level 5. Quadrants within the left half of the circle are marked for refinement, while quadrants in the right half for coarsening.

Data are mapped onto the mesh after refinement following these rules: the value of a father is inherited by its children (in case of refinement); the average sum of values of the four children is assigned to the father (in case of coarsening).

Data injection can be achieved after adaptation using an auxiliary mapper provided by an overloaded version of the adapt method. This mapper links the grid before and after adaptation.

To run: ./PABLO_example_00004
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()
{
    int iter = 0;
 
    PabloUniform pablo4(2);
 
    int idx = 0;
    pablo4.setBalance(idx,false);
 
    for (iter=1; iter<6; iter++){
        pablo4.adaptGlobalRefine();
    }
 
    double xc, yc;
    xc = yc = 0.5;
    double radius = 0.25;
 
    uint32_t nocts = pablo4.getNumOctants();
    vector<double> oct_data(nocts, 0.0);
 
    for (unsigned int i=0; i<nocts; i++){
        vector<array<double,3> > nodes = pablo4.getNodes(i);
        array<double,3> center = pablo4.getCenter(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] = (pow((center[0]-xc),2.0)+pow((center[1]-yc),2.0));
            }
        }
    }
 
    iter = 0;
    pablo4.updateConnectivity();
    pablo4.writeTest("pablo00004_iter"+to_string(static_cast<unsigned long long>(iter)), oct_data);
 
    int start = 1;
    for (iter=start; iter<start+2; iter++){
        for (unsigned int i=0; i<nocts; i++){
            vector<array<double,3> > nodes = pablo4.getNodes(i);
            array<double,3> center = pablo4.getCenter(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))){
                    if (center[0]<=xc){
 
                        pablo4.setMarker(i,1);
                    }
                    else{
 
                        pablo4.setMarker(i,-1);
                    }
                }
            }
        }
 
        vector<double> oct_data_new;
        vector<uint32_t> mapper;
        vector<bool> isghost;
        pablo4.adapt(true);
        nocts = pablo4.getNumOctants();
        oct_data_new.resize(nocts, 0.0);
 
        for (uint32_t i=0; i<nocts; i++){
            pablo4.getMapping(i, mapper, isghost);
            if (pablo4.getIsNewC(i)){
                for (int j=0; j<4; j++){
                    oct_data_new[i] += oct_data[mapper[j]]/4;
                }
            }
            else if (pablo4.getIsNewR(i)){
                oct_data_new[i] += oct_data[mapper[0]];
            }
            else{
                oct_data_new[i] += oct_data[mapper[0]];
            }
        }
 
        pablo4.updateConnectivity();
        pablo4.writeTest("pablo00004_iter"+to_string(static_cast<unsigned long long>(iter)), oct_data_new);
 
        oct_data = oct_data_new;
    }
}
 
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
}