in MultiSource/Benchmarks/DOE-ProxyApps-C++/CLAMR/mesh.cpp [2641:3610]
void Mesh::rezone_all(int icount, int jcount, vector<int> mpot, int have_state, MallocPlus &state_memory)
{
struct timeval tstart_cpu;
cpu_timer_start(&tstart_cpu);
if (! do_rezone) {
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
index.clear();
index.resize(ncells);
#ifdef _OPENMP
}
#pragma omp barrier
#endif
#ifdef _OPENMP
#pragma omp for
#endif
for (uint ic=0; ic<ncells; ic++){
index[ic]=ic;
}
#ifdef _OPENMP
#pragma omp master
#endif
cpu_timers[MESH_TIMER_REZONE_ALL] += cpu_timer_stop(tstart_cpu);
} else {
// sign for jcount is different in GPU and CPU code -- abs is a quick fix
int add_ncells = icount - abs(jcount);
#ifdef _OPENMP
#pragma omp master
#endif
cpu_counters[MESH_COUNTER_REZONE]++;
static vector<int> celltype_save;
static int new_ncells;
static int *i_old, *j_old, *level_old;
static int ifirst;
static int ilast;
static int jfirst;
static int jlast;
static int level_first;
static int level_last;
static vector<int> new_ic;
#ifdef _OPENMP
#pragma omp master
{
#endif
celltype_save.resize(ncells);
#ifdef _OPENMP
}
#pragma omp barrier
#endif
if (have_state) {
#ifdef _OPENMP
#pragma omp for
#endif
for (int ic = 0; ic < (int)ncells; ic++){
celltype_save[ic] = celltype[ic];
}
}
#ifdef _OPENMP
#pragma omp master
{
#endif
new_ncells = ncells + add_ncells;
#ifdef _OPENMP
}
#pragma omp barrier
#endif
// int ref_entry_count = 0;
if (have_state){
#ifdef _OPENMP
#pragma omp for
#endif
for (uint ic=0; ic<ncells; ic++) {
// if (mpot[ic] > 0) ref_entry_count++;
if (mpot[ic] < 0) {
// Normal cell coarsening
if (is_lower_left(i[ic],j[ic]) ) mpot[ic] = -2;
// Boundary cell case
if (celltype[ic] != REAL_CELL && is_upper_right(i[ic],j[ic]) ) mpot[ic] = -3;
}
}
}
// Initialize new variables
// int *i_old, *j_old, *level_old;
int flags = RESTART_DATA;
#ifdef HAVE_J7
if (parallel) flags = LOAD_BALANCE_MEMORY;
#endif
#ifdef _OPENMP
#pragma omp master
{
#endif
i_old = (int *)mesh_memory.memory_malloc(new_ncells, sizeof(int), "i_old", flags);
j_old = (int *)mesh_memory.memory_malloc(new_ncells, sizeof(int), "j_old", flags);
level_old = (int *)mesh_memory.memory_malloc(new_ncells, sizeof(int), "level_old", flags);
mesh_memory.memory_swap(&i, &i_old);
mesh_memory.memory_swap(&j, &j_old);
mesh_memory.memory_swap(&level, &level_old);
index.clear();
index.resize(new_ncells);
#ifdef _OPENMP
}
#pragma omp barrier
#endif
static vector<int> order; // Vector of refined mesh traversal order; set to -1 to indicate errors.
//
//vector<int> invorder(4, -1); // Vector mapping location from base index.
//int ref_entry = 0;
#ifdef _OPENMP
#pragma omp master
{
#endif
// Insert new cells into the mesh at the point of refinement.
order.resize(4, -1); // Vector of refined mesh traversal order; set to -1 to indicate errors.
ifirst = 0;
ilast = 0;
jfirst = 0;
jlast = 0;
level_first = 0;
level_last = 0;
if (parallel) {
#ifdef HAVE_MPI
MPI_Request req[12];
MPI_Status status[12];
static int prev = MPI_PROC_NULL;
static int next = MPI_PROC_NULL;
if (mype != 0) prev = mype-1;
if (mype < numpe - 1) next = mype+1;
MPI_Isend(&i_old[ncells-1], 1,MPI_INT,next,1,MPI_COMM_WORLD,req+0);
MPI_Irecv(&ifirst, 1,MPI_INT,prev,1,MPI_COMM_WORLD,req+1);
MPI_Isend(&i_old[0], 1,MPI_INT,prev,1,MPI_COMM_WORLD,req+2);
MPI_Irecv(&ilast, 1,MPI_INT,next,1,MPI_COMM_WORLD,req+3);
MPI_Isend(&j_old[ncells-1], 1,MPI_INT,next,1,MPI_COMM_WORLD,req+4);
MPI_Irecv(&jfirst, 1,MPI_INT,prev,1,MPI_COMM_WORLD,req+5);
MPI_Isend(&j_old[0], 1,MPI_INT,prev,1,MPI_COMM_WORLD,req+6);
MPI_Irecv(&jlast, 1,MPI_INT,next,1,MPI_COMM_WORLD,req+7);
MPI_Isend(&level_old[ncells-1], 1,MPI_INT,next,1,MPI_COMM_WORLD,req+8);
MPI_Irecv(&level_first, 1,MPI_INT,prev,1,MPI_COMM_WORLD,req+9);
MPI_Isend(&level_old[0], 1,MPI_INT,prev,1,MPI_COMM_WORLD,req+10);
MPI_Irecv(&level_last, 1,MPI_INT,next,1,MPI_COMM_WORLD,req+11);
MPI_Waitall(12, req, status);
#endif
}
#ifdef _OPENMP
}
#pragma omp barrier
#endif
#ifdef REZONE_NO_OPTIMIZATION
vector<int> invorder(4, -1); // Vector mapping location from base index.
for (int ic = 0, nc = 0; ic < (int)ncells; ic++)
{
if (mpot[ic] == 0 || mpot[ic] == -1000000)
{ // No change is needed; copy the old cell straight to the new mesh at this location.
index[ic] = nc;
i[nc] = i_old[ic];
j[nc] = j_old[ic];
level[nc] = level_old[ic];
nc++;
} // Complete no change needed.
else if (mpot[ic] < 0)
{ // Coarsening is needed; remove this cell and the other three and replace them with one.
index[ic] = nc;
if (mpot[ic] <= -2) {
//printf(" %d: DEBUG -- coarsening cell %d nc %d\n",mype,ic,nc);
i[nc] = i_old[ic]/2;
j[nc] = j_old[ic]/2;
level[nc] = level_old[ic] - 1;
nc++;
}
} // Coarsening complete.
else if (mpot[ic] > 0)
{ // Refinement is needed; insert four cells where once was one.
index[ic] = nc;
if (celltype[ic] == REAL_CELL)
{
set_refinement_order(&order[0], ic, ifirst, ilast, jfirst, jlast,
level_first, level_last, i_old, j_old, level_old);
// Create the cells in the correct order and orientation.
for (int ii = 0; ii < 4; ii++)
{ level[nc] = level_old[ic] + 1;
switch (order[ii])
{ case SW:
// lower left
invorder[SW] = ii;
i[nc] = i_old[ic]*2;
j[nc] = j_old[ic]*2;
nc++;
break;
case SE:
// lower right
invorder[SE] = ii;
i[nc] = i_old[ic]*2 + 1;
j[nc] = j_old[ic]*2;
nc++;
break;
case NW:
// upper left
invorder[NW] = ii;
i[nc] = i_old[ic]*2;
j[nc] = j_old[ic]*2 + 1;
nc++;
break;
case NE:
// upper right
invorder[NE] = ii;
i[nc] = i_old[ic]*2 + 1;
j[nc] = j_old[ic]*2 + 1;
nc++;
break; } } // Complete cell refinement.
} // Complete real cell refinement.
else if (celltype[ic] == LEFT_BOUNDARY) {
// lower
i[nc] = i_old[ic]*2 + 1;
j[nc] = j_old[ic]*2;
level[nc] = level_old[ic] + 1;
nc++;
// upper
i[nc] = i_old[ic]*2 + 1;
j[nc] = j_old[ic]*2 + 1;
level[nc] = level_old[ic] + 1;
nc++;
}
else if (celltype[ic] == RIGHT_BOUNDARY) {
// lower
i[nc] = i_old[ic]*2;
j[nc] = j_old[ic]*2;
level[nc] = level_old[ic] + 1;
nc++;
// upper
i[nc] = i_old[ic]*2;
j[nc] = j_old[ic]*2 + 1;
level[nc] = level_old[ic] + 1;
nc++;
}
else if (celltype[ic] == BOTTOM_BOUNDARY) {
// left
i[nc] = i_old[ic]*2;
j[nc] = j_old[ic]*2 + 1;
level[nc] = level_old[ic] + 1;
nc++;
// right
i[nc] = i_old[ic]*2 + 1;
j[nc] = j_old[ic]*2 + 1;
level[nc] = level_old[ic] + 1;
nc++;
}
else if (celltype[ic] == TOP_BOUNDARY) {
// right
i[nc] = i_old[ic]*2 + 1;
j[nc] = j_old[ic]*2;
level[nc] = level_old[ic] + 1;
nc++;
// left
i[nc] = i_old[ic]*2;
j[nc] = j_old[ic]*2;
level[nc] = level_old[ic] + 1;
nc++;
}
} // Complete refinement needed.
} // Complete addition of new cells to the mesh.
mesh_memory.memory_delete(i_old);
mesh_memory.memory_delete(j_old);
mesh_memory.memory_delete(level_old);
calc_celltype(new_ncells);
if (have_state){
flags = RESTART_DATA;
MallocPlus state_memory_old = state_memory;
malloc_plus_memory_entry *memory_item;
for (memory_item = state_memory_old.memory_entry_by_name_begin();
memory_item != state_memory_old.memory_entry_by_name_end();
memory_item = state_memory_old.memory_entry_by_name_next() ) {
//printf("DEBUG -- it.mem_name %s elsize %lu\n",memory_item->mem_name,memory_item->mem_elsize);
if (memory_item->mem_elsize == 8) {
double *state_temp_double = (double *)state_memory.memory_malloc(new_ncells, sizeof(double),
"state_temp_double", flags);
double *mem_ptr_double = (double *)memory_item->mem_ptr;
//ref_entry = 0;
for (int ic=0, nc=0; ic<(int)ncells; ic++) {
if (mpot[ic] == 0) {
state_temp_double[nc] = mem_ptr_double[ic];
nc++;
} else if (mpot[ic] < 0){
if (mpot[ic] == -2) {
int nr = nrht[ic];
int nt = ntop[ic];
int nrt = nrht[nt];
state_temp_double[nc] = (mem_ptr_double[ic] + mem_ptr_double[nr] +
mem_ptr_double[nt] + mem_ptr_double[nrt])*0.25;
nc++;
}
if (mpot[ic] == -3) {
int nl = nlft[ic];
int nb = nbot[ic];
int nlb = nlft[nb];
state_temp_double[nc] = (mem_ptr_double[ic] + mem_ptr_double[nl] +
mem_ptr_double[nb] + mem_ptr_double[nlb])*0.25;
nc++;
}
} else if (mpot[ic] > 0){
// lower left
state_temp_double[nc] = mem_ptr_double[ic];
nc++;
// lower right
state_temp_double[nc] = mem_ptr_double[ic];
nc++;
if (celltype_save[ic] == REAL_CELL){
// upper left
state_temp_double[nc] = mem_ptr_double[ic];
nc++;
// upper right
state_temp_double[nc] = mem_ptr_double[ic];
nc++;
}
}
}
state_memory.memory_replace(mem_ptr_double, state_temp_double);
} else if (memory_item->mem_elsize == 4) {
float *state_temp_float = (float *)state_memory.memory_malloc(new_ncells, sizeof(float),
"state_temp_float", flags);
float *mem_ptr_float = (float *)memory_item->mem_ptr;
for (int ic=0, nc=0; ic<(int)ncells; ic++) {
if (mpot[ic] == 0) {
state_temp_float[nc] = mem_ptr_float[ic];
nc++;
} else if (mpot[ic] < 0){
if (mpot[ic] == -2) {
int nr = nrht[ic];
int nt = ntop[ic];
int nrt = nrht[nt];
state_temp_float[nc] = (mem_ptr_float[ic] + mem_ptr_float[nr] +
mem_ptr_float[nt] + mem_ptr_float[nrt])*0.25;
nc++;
}
if (mpot[ic] == -3) {
int nl = nlft[ic];
int nb = nbot[ic];
int nlb = nlft[nb];
state_temp_float[nc] = (mem_ptr_float[ic] + mem_ptr_float[nl] +
mem_ptr_float[nb] + mem_ptr_float[nlb])*0.25;
nc++;
}
} else if (mpot[ic] > 0){
// lower left
state_temp_float[nc] = mem_ptr_float[ic];
nc++;
// lower right
state_temp_float[nc] = mem_ptr_float[ic];
nc++;
if (celltype_save[ic] == REAL_CELL){
// upper left
state_temp_float[nc] = mem_ptr_float[ic];
nc++;
// upper right
state_temp_float[nc] = mem_ptr_float[ic];
nc++;
}
}
}
state_memory.memory_replace(mem_ptr_float, state_temp_float);
}
}
}
#else
// Data parallel optimizations for thread parallel -- slows down serial
// code by about 25%
static vector<int> add_count;
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
add_count.resize(ncells);
new_ic.resize(ncells+1);
#ifdef _OPENMP
} // end master region
#pragma omp barrier
#endif
#ifdef _OPENMP
#pragma omp for
#endif
for (int ic = 0; ic < (int)ncells; ic++){
if (mpot[ic] == 0) {
add_count[ic] = 1;
} else if (mpot[ic] < 0) {
if (mpot[ic] == -2){
add_count[ic] = 1;
} else {
add_count[ic] = 0;
}
} else if (mpot[ic] > 0) {
if (celltype[ic] != REAL_CELL) {
add_count[ic] = 2;
} else {
add_count[ic] = 4;
}
}
}
#ifdef _OPENMP
#pragma omp barrier
#endif
scan (&add_count[0], &new_ic[0], ncells);
#ifdef _OPENMP
#pragma omp barrier
#endif
#ifdef _OPENMP
#pragma omp for
#endif
for (int ic = 0; ic < (int)ncells; ic++) {
vector<int> invorder(4, -1); // Vector mapping location from base index.
int nc = new_ic[ic];
if (mpot[ic] == 0)
{ // No change is needed; copy the old cell straight to the new mesh at this location.
index[ic] = nc;
i[nc] = i_old[ic];
j[nc] = j_old[ic];
level[nc] = level_old[ic];
} // Complete no change needed.
else if (mpot[ic] < 0)
{ // Coarsening is needed; remove this cell and the other three and replace them with one.
index[ic] = nc;
if (mpot[ic] <= -2) {
//printf(" %d: DEBUG -- coarsening cell %d nc %d\n",mype,ic,nc);
i[nc] = i_old[ic]/2;
j[nc] = j_old[ic]/2;
level[nc] = level_old[ic] - 1;
}
} // Coarsening complete.
else if (mpot[ic] > 0)
{ // Refinement is needed; insert four cells where once was one.
index[ic] = nc;
if (celltype[ic] == REAL_CELL)
{
int order[4];
set_refinement_order(&order[0], ic, ifirst, ilast, jfirst, jlast,
level_first, level_last, i_old, j_old, level_old);
// Create the cells in the correct order and orientation.
for (int ii = 0; ii < 4; ii++) {
level[nc] = level_old[ic] + 1;
switch (order[ii]) {
case SW:
// lower left
invorder[SW] = ii;
i[nc] = i_old[ic]*2;
j[nc] = j_old[ic]*2;
nc++;
break;
case SE:
// lower right
invorder[SE] = ii;
i[nc] = i_old[ic]*2 + 1;
j[nc] = j_old[ic]*2;
nc++;
break;
case NW:
// upper left
invorder[NW] = ii;
i[nc] = i_old[ic]*2;
j[nc] = j_old[ic]*2 + 1;
nc++;
break;
case NE:
// upper right
invorder[NE] = ii;
i[nc] = i_old[ic]*2 + 1;
j[nc] = j_old[ic]*2 + 1;
nc++;
break;
}
} // Complete cell refinement.
} // Complete real cell refinement.
else if (celltype[ic] == LEFT_BOUNDARY) {
// lower
i[nc] = i_old[ic]*2 + 1;
j[nc] = j_old[ic]*2;
level[nc] = level_old[ic] + 1;
nc++;
// upper
i[nc] = i_old[ic]*2 + 1;
j[nc] = j_old[ic]*2 + 1;
level[nc] = level_old[ic] + 1;
nc++;
}
else if (celltype[ic] == RIGHT_BOUNDARY) {
// lower
i[nc] = i_old[ic]*2;
j[nc] = j_old[ic]*2;
level[nc] = level_old[ic] + 1;
nc++;
// upper
i[nc] = i_old[ic]*2;
j[nc] = j_old[ic]*2 + 1;
level[nc] = level_old[ic] + 1;
nc++;
}
else if (celltype[ic] == BOTTOM_BOUNDARY) {
// left
i[nc] = i_old[ic]*2;
j[nc] = j_old[ic]*2 + 1;
level[nc] = level_old[ic] + 1;
nc++;
// right
i[nc] = i_old[ic]*2 + 1;
j[nc] = j_old[ic]*2 + 1;
level[nc] = level_old[ic] + 1;
nc++;
}
else if (celltype[ic] == TOP_BOUNDARY) {
// right
i[nc] = i_old[ic]*2 + 1;
j[nc] = j_old[ic]*2;
level[nc] = level_old[ic] + 1;
nc++;
// left
i[nc] = i_old[ic]*2;
j[nc] = j_old[ic]*2;
level[nc] = level_old[ic] + 1;
nc++;
}
} // Complete refinement needed.
} // Complete addition of new cells to the mesh.
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
mesh_memory.memory_delete(i_old);
mesh_memory.memory_delete(j_old);
mesh_memory.memory_delete(level_old);
#ifdef _OPENMP
} // end master region
#endif
calc_celltype_threaded(new_ncells);
if (have_state){
static MallocPlus state_memory_old;
static malloc_plus_memory_entry *memory_begin;
static malloc_plus_memory_entry *memory_end;
static malloc_plus_memory_entry *memory_next;
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
state_memory_old = state_memory;
memory_begin = state_memory_old.memory_entry_by_name_begin();
memory_end = state_memory_old.memory_entry_by_name_end();
#ifdef _OPENMP
} // end master region
#pragma omp barrier
#endif
for (malloc_plus_memory_entry *memory_item = memory_begin;
memory_item != memory_end;
memory_item = memory_next ) {
//ref_entry = 0;
//printf("DEBUG -- memory_item->mem_name %s elsize %lu\n",memory_item->mem_name,memory_item->mem_elsize);
if (memory_item->mem_elsize == 8) {
static double *state_temp_double, *mem_ptr_double;
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
state_temp_double = (double *)state_memory.memory_malloc(new_ncells, sizeof(double),
"state_temp_double", flags);
mem_ptr_double = (double *)memory_item->mem_ptr;
#ifdef _OPENMP
} // end master region
#pragma omp barrier
#endif
//ref_entry = 0;
#ifdef _OPENMP
#pragma omp for
#endif
for (int ic=0; ic<(int)ncells; ic++) {
int nc = new_ic[ic];
if (mpot[ic] == 0) {
state_temp_double[nc] = mem_ptr_double[ic];
} else if (mpot[ic] < 0){
if (mpot[ic] == -2) {
int nr = nrht[ic];
int nt = ntop[ic];
int nrt = nrht[nt];
state_temp_double[nc] = (mem_ptr_double[ic] + mem_ptr_double[nr] +
mem_ptr_double[nt] + mem_ptr_double[nrt])*0.25;
}
if (mpot[ic] == -3) {
int nl = nlft[ic];
int nb = nbot[ic];
int nlb = nlft[nb];
state_temp_double[nc] = (mem_ptr_double[ic] + mem_ptr_double[nl] +
mem_ptr_double[nb] + mem_ptr_double[nlb])*0.25;
}
} else if (mpot[ic] > 0){
// lower left
state_temp_double[nc] = mem_ptr_double[ic];
nc++;
// lower right
state_temp_double[nc] = mem_ptr_double[ic];
nc++;
if (celltype_save[ic] == REAL_CELL){
// upper left
state_temp_double[nc] = mem_ptr_double[ic];
nc++;
// upper right
state_temp_double[nc] = mem_ptr_double[ic];
nc++;
}
}
} // end cell loop
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
state_memory.memory_replace(mem_ptr_double, state_temp_double);
#ifdef _OPENMP
} // end master region
#pragma omp barrier
#endif
} else if (memory_item->mem_elsize == 4) {
static float *state_temp_float, *mem_ptr_float;
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
state_temp_float = (float *)state_memory.memory_malloc(new_ncells, sizeof(float),
"state_temp_float", flags);
mem_ptr_float = (float *)memory_item->mem_ptr;
#ifdef _OPENMP
} // end master region
#pragma omp barrier
#endif
#ifdef _OPENMP
#pragma omp for
#endif
for (int ic=0; ic<(int)ncells; ic++) {
int nc = new_ic[ic];
if (mpot[ic] == 0) {
state_temp_float[nc] = mem_ptr_float[ic];
} else if (mpot[ic] < 0){
if (mpot[ic] == -2) {
int nr = nrht[ic];
int nt = ntop[ic];
int nrt = nrht[nt];
state_temp_float[nc] = (mem_ptr_float[ic] + mem_ptr_float[nr] +
mem_ptr_float[nt] + mem_ptr_float[nrt])*0.25;
}
if (mpot[ic] == -3) {
int nl = nlft[ic];
int nb = nbot[ic];
int nlb = nlft[nb];
state_temp_float[nc] = (mem_ptr_float[ic] + mem_ptr_float[nl] +
mem_ptr_float[nb] + mem_ptr_float[nlb])*0.25;
}
} else if (mpot[ic] > 0){
// lower left
state_temp_float[nc] = mem_ptr_float[ic];
nc++;
// lower right
state_temp_float[nc] = mem_ptr_float[ic];
nc++;
if (celltype_save[ic] == REAL_CELL){
// upper left
state_temp_float[nc] = mem_ptr_float[ic];
nc++;
// upper right
state_temp_float[nc] = mem_ptr_float[ic];
nc++;
}
}
} // end cell loop
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
state_memory.memory_replace(mem_ptr_float, state_temp_float);
#ifdef _OPENMP
} // end master region
#pragma omp barrier
#endif
} // mem elem size 4 bytes
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
memory_next = state_memory_old.memory_entry_by_name_next();
#ifdef _OPENMP
} // end master region
#pragma omp barrier
#endif
} // memory item iteration
} // if have state
// End of data parallel optimizations
#endif
if (neighbor_remap) {
int flags = 0;
static int *nlft_old, *nrht_old, *nbot_old, *ntop_old;
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
nlft_old = (int *)mesh_memory.memory_malloc(new_ncells, sizeof(int), "nlft_old", flags);
nrht_old = (int *)mesh_memory.memory_malloc(new_ncells, sizeof(int), "nrht_old", flags);
nbot_old = (int *)mesh_memory.memory_malloc(new_ncells, sizeof(int), "nbot_old", flags);
ntop_old = (int *)mesh_memory.memory_malloc(new_ncells, sizeof(int), "ntop_old", flags);
#ifdef _OPENMP
} // end master region
#pragma omp barrier
#endif
flags = RESTART_DATA;
#ifdef _OPENMP
#pragma omp for
#endif
for (int ic = 0; ic < new_ncells; ic++){
nlft_old[ic] = -1;
nrht_old[ic] = -1;
nbot_old[ic] = -1;
ntop_old[ic] = -1;
}
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
mesh_memory.memory_swap(&nlft, &nlft_old);
mesh_memory.memory_swap(&nrht, &nrht_old);
mesh_memory.memory_swap(&nbot, &nbot_old);
mesh_memory.memory_swap(&ntop, &ntop_old);
#ifdef _OPENMP
} // end master region
#pragma omp barrier
#endif
#ifdef _OPENMP
#pragma omp for
#endif
for (int ic = 0; ic < (int)ncells; ic++){
int nc = index[ic];
if (mpot[ic] == 0){
if (nlft_old[ic] < (int)ncells && nlft_old[ic] >= 0){
nlft[nc] = (mpot[nlft_old[ic]] == 0) ? index[nlft_old[ic]] : -1;
}
if (nrht_old[ic] < (int)ncells && nrht_old[ic] >= 0){
nrht[nc] = (mpot[nrht_old[ic]] == 0) ? index[nrht_old[ic]] : -1;
}
if (nbot_old[ic] < (int)ncells && nbot_old[ic] >= 0){
nbot[nc] = (mpot[nbot_old[ic]] == 0) ? index[nbot_old[ic]] : -1;
}
if (ntop_old[ic] < (int)ncells && ntop_old[ic] >= 0){
ntop[nc] = (mpot[ntop_old[ic]] == 0) ? index[ntop_old[ic]] : -1;
}
} else if (mpot[ic] <= -2) {
nlft[nc] = -1;
nrht[nc] = -1;
nbot[nc] = -1;
ntop[nc] = -1;
} else if (mpot[ic] > 0){
nlft[nc] = -1;
nlft[nc+1] = -1;
nrht[nc] = -1;
nrht[nc+1] = -1;
nbot[nc] = -1;
nbot[nc+1] = -1;
ntop[nc] = -1;
ntop[nc+1] = -1;
if (celltype[nc] == REAL_CELL){
nlft[nc+2] = -1;
nlft[nc+3] = -1;
nrht[nc+2] = -1;
nrht[nc+3] = -1;
nbot[nc+2] = -1;
nbot[nc+3] = -1;
ntop[nc+2] = -1;
ntop[nc+3] = -1;
}
}
if (mpot[ic] > 0){
nc++;
switch(celltype[nc]){
case LEFT_BOUNDARY:
nlft[nc] = nc;
break;
case RIGHT_BOUNDARY:
nrht[nc] = nc;
break;
case BOTTOM_BOUNDARY:
nbot[nc] = nc;
break;
case TOP_BOUNDARY:
ntop[nc] = nc;
break;
}
}
}
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
nlft_old = (int *)mesh_memory.memory_delete(nlft_old);
nrht_old = (int *)mesh_memory.memory_delete(nrht_old);
nbot_old = (int *)mesh_memory.memory_delete(nbot_old);
ntop_old = (int *)mesh_memory.memory_delete(ntop_old);
#ifdef _OPENMP
} // end master region
#pragma omp barrier
#endif
} else {
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
nlft = (int *)mesh_memory.memory_delete(nlft);
nrht = (int *)mesh_memory.memory_delete(nrht);
nbot = (int *)mesh_memory.memory_delete(nbot);
ntop = (int *)mesh_memory.memory_delete(ntop);
#ifdef _OPENMP
} // end master region
#pragma omp barrier
#endif
}
#ifdef _OPENMP
#pragma omp barrier
#pragma omp master
{
#endif
//ncells = nc;
#ifdef HAVE_MPI
if (parallel) {
MPI_Allgather(&new_ncells, 1, MPI_INT, &nsizes[0], 1, MPI_INT, MPI_COMM_WORLD);
ndispl[0]=0;
for (int ip=1; ip<numpe; ip++){
ndispl[ip] = ndispl[ip-1] + nsizes[ip-1];
}
noffset=ndispl[mype];
ncells_global = ndispl[numpe-1]+nsizes[numpe-1];
}
#endif
cpu_timers[MESH_TIMER_REZONE_ALL] += cpu_timer_stop(tstart_cpu);
#ifdef _OPENMP
} // end master region
#pragma omp barrier
#endif
} // if do_rezone
}