vtkStructuredAMRGridConnectivity.cxx 75.6 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
/*=========================================================================

 Program:   Visualization Toolkit
 Module:    vtkStructuredAMRGridConnectivity.cxx

 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
 All rights reserved.
 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

 This software is distributed WITHOUT ANY WARRANTY; without even
 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 PURPOSE.  See the above copyright notice for more information.

 =========================================================================*/
#include "vtkStructuredAMRGridConnectivity.h"

// VTK Includes
#include "vtkObjectFactory.h"
#include "vtkStructuredData.h"
#include "vtkStructuredGridConnectivity.h"
#include "vtkStructuredNeighbor.h"

// C++ includes
#include <cmath>
#include <cstdlib>

//------------------------------------------------------------------------------
#define IMIN(ext) ext[0]
#define IMAX(ext) ext[1]
#define JMIN(ext) ext[2]
#define JMAX(ext) ext[3]
#define KMIN(ext) ext[4]
#define KMAX(ext) ext[5]

namespace AMRBlockFace {
  enum
37
  {
38 39 40 41 42 43 44
    FRONT         = 0, // (+k diretion)
    BACK          = 1, // (-k direction)
    RIGHT         = 2, // (+i direction)
    LEFT          = 3, // (-i direction)
    TOP           = 4, // (+j direction)
    BOTTOM        = 5, // (-j direction)
    NOT_ON_BLOCK_FACE = 6
45
  };
46 47 48 49 50 51 52 53 54 55 56
}

vtkStandardNewMacro( vtkStructuredAMRGridConnectivity );

//-----------------------------------------------------------------------------
vtkStructuredAMRGridConnectivity::vtkStructuredAMRGridConnectivity()
{
  this->DataDimension      = 0;
  this->DataDescription    = VTK_EMPTY;
  this->NumberOfGrids      = 0;
  this->MaxLevel           = -1;
57
  this->RefinementRatio    = -1;
58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89
  this->NumberOfLevels     = 0;
  this->BalancedRefinement = true;
  this->CellCentered       = true;
  this->NodeCentered       = false;

  IMIN(this->WholeExtent) =
  JMIN(this->WholeExtent) =
  KMIN(this->WholeExtent) = VTK_INT_MAX;

  IMAX(this->WholeExtent) =
  JMAX(this->WholeExtent) =
  KMAX(this->WholeExtent) = VTK_INT_MIN;
}

//-----------------------------------------------------------------------------
vtkStructuredAMRGridConnectivity::~vtkStructuredAMRGridConnectivity()
{
  this->AMRHierarchy.clear();
  this->GridExtents.clear();
  this->GridLevels.clear();
  this->Neighbors.clear();
}

//-----------------------------------------------------------------------------
void vtkStructuredAMRGridConnectivity::PrintSelf(
    std::ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf( os, indent );
  os << "=====================\n";
  os << "DATA DIMENSION: " << this->DataDimension << std::endl;
  os << "WHOLE EXTENT: [";
  for( int i=0; i < 6; ++i )
90
  {
91
    os << this->WholeExtent[i] << " ";
92
  }
93 94 95 96
  os << "]\n";
  os << "TOTAL NUMBER OF LEVELS:" << this->NumberOfLevels << std::endl;
  os << "TOTAL NUMBER OF GRIDS:"  << this->NumberOfGrids  << std::endl;
  if( this->HasConstantRefinementRatio() )
97
  {
98
    os << "CONSTANT REFINEMENT RATIO: " << this->RefinementRatio << std::endl;
99
  }
100
  else
101
  {
102
    os << "VARIABLE REFINEMENT RATIO\n";
103
  }
104 105 106 107

  int gridExtent[6];
  int neiExtent[6];
  for( unsigned int gridId=0; gridId < this->NumberOfGrids; ++gridId )
108
  {
109 110 111 112 113 114 115
    os << "=====================\n";
    os << "GRID["  << gridId << "] ";
    os << "LEVEL=" << this->GetGridLevel( gridId ) << " ";
    os << "EXTENT: ";
    this->GetGridExtent(gridId,gridExtent);
    this->PrintExtent(os,gridExtent);
    os << std::endl;
116
    if( !this->GhostedExtents.empty() )
117
    {
118 119 120 121 122 123 124
      assert("pre: ghosted extents vector is not properly allocated" &&
             (this->GhostedExtents.size()/6 == this->NumberOfGrids) );
      os << "GHOSTED EXTENT: ";
      int ghostedExt[6];
      this->GetGhostedExtent( gridId, ghostedExt );
      this->PrintExtent(os,ghostedExt);
      os << std::endl;
125
    }
126 127 128 129 130 131 132

    os << std::endl;
    os << "Connecting faces: "
       << this->GetNumberOfConnectingBlockFaces( gridId ) << " ";

    os << "[ ";
    if( this->HasBlockConnection( gridId, AMRBlockFace::FRONT ) )
133
    {
134
      os << "FRONT(+k) ";
135
    }
136
    if( this->HasBlockConnection( gridId, AMRBlockFace::BACK ) )
137
    {
138
      os << "BACK(-k) ";
139
    }
140
    if( this->HasBlockConnection( gridId, AMRBlockFace::RIGHT ) )
141
    {
142
      os << "RIGHT(+i) ";
143
    }
144
    if( this->HasBlockConnection( gridId, AMRBlockFace::LEFT ) )
145
    {
146
      os << "LEFT(-i) ";
147
    }
148
    if( this->HasBlockConnection( gridId, AMRBlockFace::TOP) )
149
    {
150
      os << "TOP(+j) ";
151
    }
152
    if( this->HasBlockConnection( gridId, AMRBlockFace::BOTTOM) )
153
    {
154
      os << "BOTTOM(-j) ";
155
    }
156 157 158 159 160 161 162
    os << "] ";
    os << std::endl;

    os << "NUMBER OF NEIGHBORS: " << this->Neighbors[gridId].size();
    os << std::endl << std::endl;

    for( unsigned int nei=0; nei < this->Neighbors[gridId].size(); ++nei)
163
    {
164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195
      os << "\t-----------------------------" << std::endl;
      os << "\tNEIGHBOR[" << nei << "] ";
      os << "ID=" << this->Neighbors[gridId][nei].NeighborID << " ";
      os << "LEVEL=" << this->Neighbors[gridId][nei].NeighborLevel << " ";
      os << "EXTENT=";
      this->GetGridExtent(this->Neighbors[gridId][nei].NeighborID,neiExtent);
      this->PrintExtent(os,neiExtent);
      os << " RELATIONSHIP=";
      os << this->Neighbors[gridId][nei].GetRelationShipString();
      os << std::endl;

      os << "\tGRID OVERLAP EXTENT=";
      this->PrintExtent(os,this->Neighbors[gridId][nei].GridOverlapExtent);
      os << "NEI OVERLAP EXTENT=";
      this->PrintExtent(os, this->Neighbors[gridId][nei].OverlapExtent);
      os << std::endl;

      os << "\tORIENTATION: (";
      os << this->Neighbors[gridId][nei].Orientation[0];
      os << ", ";
      os << this->Neighbors[gridId][nei].Orientation[1];
      os << ", ";
      os << this->Neighbors[gridId][nei].Orientation[2];
      os << ")\n";
      os << std::endl << std::endl;

      os << "\tRCVEXTENT=";
      this->PrintExtent(os,this->Neighbors[gridId][nei].RcvExtent);
      os << std::endl;
      os << "\tSNDEXTENT=";
      this->PrintExtent(os,this->Neighbors[gridId][nei].SendExtent);
      os << std::endl << std::endl;
196
    } // END for all neighbors
197

198
  } // END for all grids
199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214
}

//-----------------------------------------------------------------------------
void vtkStructuredAMRGridConnectivity::Initialize(
    const unsigned int NumLevels,
    const unsigned int N, const int rr)
{
  this->NumberOfLevels = NumLevels;
  this->RefinementRatio = rr;
  this->SetNumberOfGrids( N );
}

//-----------------------------------------------------------------------------
void vtkStructuredAMRGridConnectivity::SetNumberOfGrids(
    const unsigned int N )
{
Andrew Maclean's avatar
Andrew Maclean committed
215
  if (N == 0)
216
  {
Andrew Maclean's avatar
Andrew Maclean committed
217 218
    vtkErrorMacro("Number of grids cannot be 0.");
    return;
219
  }
220 221 222 223 224 225 226 227 228
  this->NumberOfGrids = N;
  this->AllocateUserRegisterDataStructures();

  this->GridExtents.resize( 6*N );
  this->GridLevels.resize( N );
  this->Neighbors.resize( N );
  this->BlockTopology.resize( N );

  if( !this->HasConstantRefinementRatio() )
229
  {
230
    this->RefinementRatios.resize( this->NumberOfLevels,-1 );
231
  }
232 233 234 235 236 237 238 239 240 241 242 243 244 245 246
}

//-----------------------------------------------------------------------------
void vtkStructuredAMRGridConnectivity::SetBlockTopology(
    const int gridID)
{
  assert("pre: gridID is out-of-bounds!" &&
          (gridID >= 0) &&
          (gridID < static_cast<int>(this->NumberOfGrids)));

  int gridExtent[6];
  this->GetCoarsenedExtent(gridID,this->GridLevels[gridID],0,gridExtent);

  // Check in IMIN
  if( gridExtent[0] > this->WholeExtent[0] )
247
  {
248
    this->AddBlockConnection( gridID, AMRBlockFace::LEFT);
249
  }
250 251 252

  // Check in IMAX
  if( gridExtent[1] < this->WholeExtent[1] )
253
  {
254
    this->AddBlockConnection( gridID, AMRBlockFace::RIGHT);
255
  }
256 257 258

  // Check in JMIN
  if( gridExtent[2] > this->WholeExtent[2] )
259
  {
260
    this->AddBlockConnection( gridID, AMRBlockFace::BOTTOM );
261
  }
262 263 264

  // Check in JMAX
  if( gridExtent[3] < this->WholeExtent[3] )
265
  {
266
    this->AddBlockConnection( gridID, AMRBlockFace::TOP );
267
  }
268 269 270

  // Check in KMIN
  if( gridExtent[4] > this->WholeExtent[4] )
271
  {
272
    this->AddBlockConnection( gridID, AMRBlockFace::BACK );
273
  }
274 275 276

  // Check in KMAX
  if( gridExtent[5] < this->WholeExtent[5] )
277
  {
278
    this->AddBlockConnection( gridID, AMRBlockFace::FRONT );
279
  }
280 281 282 283 284 285 286 287 288 289 290
}

//-----------------------------------------------------------------------------
void vtkStructuredAMRGridConnectivity::ComputeNeighbors()
{
  // STEP 0: Compute the whole extent w.r.t. level 0 which also computes the
  // data-description and dimension of the data.
  this->ComputeWholeExtent();

  // STEP 1: Establish neighbor relation between grids in the AMR hierarchy
  for( unsigned int i=0; i < this->NumberOfGrids; ++i )
291
  {
292 293 294
    this->SetBlockTopology( i );

    for( unsigned int j=i+1; j < this->NumberOfGrids; ++j )
295
    {
296
      this->EstablishNeighbors(i,j);
297
    } // END for all j
298 299 300

    this->FillGhostArrays(
        i,this->GridPointGhostArrays[i],this->GridCellGhostArrays[i]);
301
  } // END for all i
302 303 304 305 306 307 308 309 310 311 312 313 314 315
}

//-----------------------------------------------------------------------------
void vtkStructuredAMRGridConnectivity::GetGhostedExtent(
        const int gridID, int ext[6])
{
  assert("pre: grid ID is out-of-bounds!" &&
       (gridID >= 0) && (gridID < static_cast<int>(this->NumberOfGrids)));
  assert("pre: ghosted-extents vector has not been allocated" &&
       (this->NumberOfGrids==this->GhostedExtents.size()/6));
  assert("pre: Number of ghost layers should not be 0" &&
       (this->NumberOfGhostLayers > 0) );

  for( int i=0; i < 6; ++i )
316
  {
317
    ext[ i ] = this->GhostedExtents[gridID*6+i];
318
  }
319 320 321 322 323 324 325 326 327 328 329 330 331 332
}

//-----------------------------------------------------------------------------
void vtkStructuredAMRGridConnectivity::SetGhostedExtent(
        const int gridID, int ext[6])
{
  assert("pre: grid ID is out-of-bounds!" &&
       (gridID >= 0) && (gridID < static_cast<int>(this->NumberOfGrids)));
  assert("pre: ghosted-extents vector has not been allocated" &&
       (this->NumberOfGrids==this->GhostedExtents.size()/6));
  assert("pre: Number of ghost layers should not be 0" &&
       (this->NumberOfGhostLayers > 0) );

  for( int i=0; i < 6; ++i )
333
  {
334
    this->GhostedExtents[gridID*6+i] = ext[ i ] ;
335
  }
336 337 338 339 340 341
}

//-----------------------------------------------------------------------------
void vtkStructuredAMRGridConnectivity::CreateGhostLayers(const int N)
{
  if( N==0 )
342
  {
343 344 345
    vtkWarningMacro(
        "N=0 ghost layers requested! No ghost layers will be created");
    return;
346
  }
347 348 349 350 351 352

  this->NumberOfGhostLayers += N;
  this->AllocateInternalDataStructures();
  this->GhostedExtents.resize( 6*this->NumberOfGrids );

  for(unsigned int i=0; i < this->NumberOfGrids; ++i)
353
  {
354 355 356 357 358 359
    this->CreateGhostedExtent( i, N );
    this->CreateGhostedMaskArrays( i );
    this->ComputeNeighborSendAndRcvExtent(i,N);
    this->InitializeGhostData( i );
    this->TransferRegisteredDataToGhostedData( i );
    this->TransferGhostDataFromNeighbors( i );
360
  } // END for all grids
361 362 363 364 365 366 367 368 369
}

//-----------------------------------------------------------------------------
void vtkStructuredAMRGridConnectivity::InitializeGhostData(
    const int gridID)
{
  assert( "pre: gridID is out-of-bounds!" &&
          (gridID >= 0) && (gridID < static_cast<int>(this->NumberOfGrids)));
  assert( "pre: Grid has no registered point data!" &&
370
          (this->GridPointData[gridID] != nullptr) );
371
  assert( "pre: Grid has no registered cell data!" &&
372
          (this->GridCellData[gridID] != nullptr) );
373 374 375 376 377 378 379

  // STEP 0: Get the ghosted grid extent
  int ghostedExtent[6];
  this->GetGhostedExtent( gridID, ghostedExtent );

  // STEP 1: Get the number of nodes/cells in the ghosted extent
  int numNodes =
380
    vtkStructuredData::GetNumberOfPoints(ghostedExtent,this->DataDescription);
381 382 383 384 385 386 387 388
  int numCells =
    vtkStructuredData::GetNumberOfCells(ghostedExtent,this->DataDescription);

  // NOTE: For AMR we currently only support uniform AMR, so there is no need
  // to allocate the GhostedGridPoints

  // STEP 2: Allocate point data, if node-centered is true
  if( this->GetNodeCentered() )
389
  {
390 391 392 393 394 395
    assert( "pre: GhostedPointData vector has not been properly allocated!" &&
         (this->NumberOfGrids==this->GhostedGridPointData.size() ) );

    this->GhostedGridPointData[ gridID ] = vtkPointData::New();
    vtkPointData *PD = this->GridPointData[gridID];
    for(int array=0; array < PD->GetNumberOfArrays(); ++array)
396
    {
397 398
      int dataType = PD->GetArray( array )->GetDataType();
      vtkDataArray *dataArray = vtkDataArray::CreateDataArray( dataType );
399
      assert( "Cannot create data array" && (dataArray != nullptr) );
400 401 402 403 404 405 406 407

      dataArray->SetName(PD->GetArray(array)->GetName());
      dataArray->SetNumberOfComponents(
          PD->GetArray(array)->GetNumberOfComponents());
      dataArray->SetNumberOfTuples( numNodes );

      this->GhostedGridPointData[ gridID ]->AddArray( dataArray );
      dataArray->Delete();
408 409
    } // END for all node arrays
  } // END if node-centered data-set
410 411 412

  // STEP 3: Allocate cell data
  if( this->GetCellCentered() )
413
  {
414 415 416 417 418
    assert("pre: GhostedCellData vector has not been properly allocated!" &&
           (this->NumberOfGrids==this->GhostedGridCellData.size() ) );
    this->GhostedGridCellData[ gridID ] = vtkCellData::New();
    vtkCellData *CD = this->GridCellData[gridID];
    for(int array=0; array < CD->GetNumberOfArrays(); ++array)
419
    {
420 421
      int dataType = CD->GetArray( array )->GetDataType();
      vtkDataArray *dataArray = vtkDataArray::CreateDataArray( dataType );
422
      assert( "Cannot create data array" && (dataArray != nullptr) );
423 424 425 426 427 428 429 430

      dataArray->SetName(CD->GetArray(array)->GetName());
      dataArray->SetNumberOfComponents(
          CD->GetArray(array)->GetNumberOfComponents());
      dataArray->SetNumberOfTuples(numCells);

      this->GhostedGridCellData[ gridID ]->AddArray(dataArray);
      dataArray->Delete();
431 432
    } // END for all cell arrays
  } // END if cell-centered data-set
433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456
}

//-----------------------------------------------------------------------------
void
vtkStructuredAMRGridConnectivity::TransferRegisteredDataToGhostedData(
    const int gridID)
{
  assert("pre: grid ID is out-of-bounds!" &&
       (gridID >= 0) && (gridID < static_cast<int>(this->NumberOfGrids)));

  // NOTE: For AMR we only support uniform grids, so we only transfer fields,
  // i.e., PointData and CellData here.

  // STEP 0: Get the registered grid extent
  int registeredExtent[6];
  this->GetGridExtent(gridID,registeredExtent);

  // STEP 1: Get the ghosted grid extent
  int ghostedExtent[6];
  this->GetGhostedExtent(gridID,ghostedExtent);

  // STEP 2: Get corresponding registered and ghosted cell extents
  int registeredCellExtent[6];
  int ghostedCellExtent[6];
457
  vtkStructuredData::GetCellExtentFromPointExtent(
458
      registeredExtent,registeredCellExtent,this->DataDescription);
459
  vtkStructuredData::GetCellExtentFromPointExtent(
460 461 462 463 464
      ghostedExtent,ghostedCellExtent,this->DataDescription);

  // STEP 3: Loop over registered grid extent
  int ijk[3];
  for(int i=IMIN(registeredExtent); i <= IMAX(registeredExtent); ++i)
465
  {
466
    for(int j=JMIN(registeredExtent); j <= JMAX(registeredExtent); ++j)
467
    {
468
      for(int k=KMIN(registeredExtent); k <= KMAX(registeredExtent); ++k)
469
      {
470 471 472
        ijk[0]=i; ijk[1]=j; ijk[2]=k;

        if( this->GetNodeCentered() )
473
        {
474 475 476 477 478 479 480 481 482 483 484 485 486
          // Compute the source index to the registered data
          vtkIdType sourcePntIdx =
              vtkStructuredData::ComputePointIdForExtent(
                  registeredExtent,ijk,this->DataDescription);

          // Compute the target index to the ghosted data
          vtkIdType targetPntIdx =
              vtkStructuredData::ComputePointIdForExtent(
                  ghostedExtent,ijk,this->DataDescription);

          this->CopyFieldData(
              this->GridPointData[gridID],sourcePntIdx,
              this->GhostedGridPointData[gridID],targetPntIdx);
487
        }
488 489

        if( this->IsNodeWithinExtent(i,j,k,registeredCellExtent) )
490
        {
491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508
          // Compute the source cell idx. Note, since we are passing to
          // ComputePointIdForExtent a cell extent, this is a cell id, not
          // a point id.
          vtkIdType sourceCellIdx =
              vtkStructuredData::ComputePointIdForExtent(
                  registeredCellExtent, ijk, this->DataDescription );

          // Compute the target cell idx. Note, since we are passing to
          // ComputePointIdForExtent a cell extent, this is a cell id, not
          // a point id.
          vtkIdType targetCellIdx =
              vtkStructuredData::ComputePointIdForExtent(
                  ghostedCellExtent, ijk, this->DataDescription );

          // Transfer cell data from the registered grid to the ghosted grid
          this->CopyFieldData(
              this->GridCellData[gridID], sourceCellIdx,
              this->GhostedGridCellData[gridID], targetCellIdx );
509
        } // END if within the cell extent
510

511 512 513
      } // END for all k
    } // END for all j
  } // END for all i
514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533
}

//-----------------------------------------------------------------------------
void
vtkStructuredAMRGridConnectivity::TransferLocalNodeCenteredNeighborData(
    const int vtkNotUsed(gridID),
    vtkStructuredAMRNeighbor& vtkNotUsed(nei))
{
  vtkErrorMacro("Node-centered AMR datasets are currently not supported!");
}

//-----------------------------------------------------------------------------
void
vtkStructuredAMRGridConnectivity::GetLocalCellCentersAtSameLevel(
      const int gridID, vtkStructuredAMRNeighbor &nei)
{
  // STEP 0: Get the grid's extent and cell extent
  int RegisteredGridExtent[6];
  this->GetGridExtent( gridID, RegisteredGridExtent );
  int RegisteredGridCellExtent[6];
534
  vtkStructuredData::GetCellExtentFromPointExtent(
535 536 537 538 539 540
      RegisteredGridExtent,RegisteredGridCellExtent,this->DataDescription);

  // STEP 1: Get the grid's ghosted extent and cell extent
  int GhostedGridExtent[6];
  this->GetGhostedExtent( gridID, GhostedGridExtent );
  int GhostedCellExtent[6];
541
  vtkStructuredData::GetCellExtentFromPointExtent(
542 543 544 545 546 547 548
      GhostedGridExtent,GhostedCellExtent,this->DataDescription);


  // STEP 2: Get the neighbor's extent and cell extent
  int NeighborExtent[6];
  this->GetGridExtent( nei.NeighborID, NeighborExtent );
  int NeighborCellExtent[6];
549
  vtkStructuredData::GetCellExtentFromPointExtent(
550 551 552 553
      NeighborExtent,NeighborCellExtent,this->DataDescription);

  // STEP 3: Get RcvCell extent
  int RcvCellExtent[6];
554
  vtkStructuredData::GetCellExtentFromPointExtent(
555 556 557 558 559 560
      const_cast<int*>(nei.RcvExtent),RcvCellExtent);

  // STEP 4: Loop through the RcvCellExtent and copy values iff a higher res
  // value does not exist.
  int ijk[3];
  for (int i=IMIN(RcvCellExtent); i <= IMAX(RcvCellExtent); ++i)
561
  {
562
    for (int j=JMIN(RcvCellExtent); j <= JMAX(RcvCellExtent); ++j)
563
    {
564
      for (int k=KMIN(RcvCellExtent); k <= KMAX(RcvCellExtent); ++k)
565
      {
566 567 568 569
        ijk[0]=i; ijk[1]=j; ijk[2]=k;

        if( this->IsNodeWithinExtent(i,j,k, NeighborCellExtent)
            && !this->IsNodeWithinExtent(i,j,k, RegisteredGridCellExtent))
570
        {
571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586
          // Sanity check!
          assert("pre: RcvExtent is outside the GhostExtent!" &&
              this->IsNodeWithinExtent(i,j,k,GhostedGridExtent));
          assert("pre: RcvExtent is outside the NeighborExtent" &&
              this->IsNodeWithinExtent(i,j,k,NeighborExtent));

          // Compute the source & target index
          // Note: Since these indices are computed from a cell extent they
          // correspond to a cell index.
          vtkIdType sourceIdx = vtkStructuredData::ComputePointIdForExtent(
              NeighborCellExtent, ijk, this->DataDescription);

          vtkIdType targetIdx = vtkStructuredData::ComputePointIdForExtent(
              GhostedCellExtent, ijk, this->DataDescription);

          if(this->CellCenteredDonorLevel[gridID][targetIdx]< nei.NeighborLevel)
587
          {
588 589 590 591
            this->CopyFieldData(
                this->GridCellData[nei.NeighborID],sourceIdx,
                this->GhostedGridCellData[gridID],targetIdx);
            this->CellCenteredDonorLevel[gridID][targetIdx]=nei.NeighborLevel;
592 593
          } // END if this is a finer solution
        } // END if this is a ghost cell
594

595 596 597
      } // END for all k
    } // END for all j
  } // END for all i
598 599 600 601 602 603 604 605 606 607 608 609 610 611 612

}

//-----------------------------------------------------------------------------
void
vtkStructuredAMRGridConnectivity::GetLocalCellCentersFromCoarserLevel(
    const int gridID, vtkStructuredAMRNeighbor &nei)
{
  assert( "pre: Expected a coarser neighbor" &&
          (nei.NeighborLevel < nei.GridLevel) );

  // STEP 0: Get the grid's extent and cell extent
  int RegisteredGridExtent[6];
  this->GetGridExtent( gridID, RegisteredGridExtent );
  int RegisteredGridCellExtent[6];
613
  vtkStructuredData::GetCellExtentFromPointExtent(
614 615 616 617 618 619
      RegisteredGridExtent,RegisteredGridCellExtent,this->DataDescription);

  // STEP 1: Get the grid's ghosted extent and cell extent
  int GhostedGridExtent[6];
  this->GetGhostedExtent( gridID, GhostedGridExtent );
  int GhostedCellExtent[6];
620
  vtkStructuredData::GetCellExtentFromPointExtent(
621 622 623 624 625 626 627
      GhostedGridExtent,GhostedCellExtent,this->DataDescription);


  // STEP 3: Get the neighbor's extent and cell extent
  int NeighborExtent[6];
  this->GetGridExtent( nei.NeighborID, NeighborExtent );
  int NeighborCellExtent[6];
628
  vtkStructuredData::GetCellExtentFromPointExtent(
629 630 631 632 633 634
      NeighborExtent,NeighborCellExtent,this->DataDescription);

  // STEP 4: Get RcvCell extent
  int rcvDataDescription = vtkStructuredData::GetDataDescriptionFromExtent(
                              const_cast<int*>(nei.RcvExtent));
  int RcvCellExtent[6];
635
  vtkStructuredData::GetCellExtentFromPointExtent(
636 637 638 639 640
      const_cast<int*>(nei.RcvExtent),RcvCellExtent);

  // STEP 5: Loop through the rcv cell extent and fill ghost regions
  int ijk[3];
  for(int i=IMIN(RcvCellExtent); i <= IMAX(RcvCellExtent); ++i)
641
  {
642
    for(int j=JMIN(RcvCellExtent); j <= JMAX(RcvCellExtent); ++j)
643
    {
644
      for(int k=KMIN(RcvCellExtent); k <= KMAX(RcvCellExtent); ++k)
645
      {
646 647 648 649 650 651 652 653 654 655 656 657 658 659
        ijk[0]=i; ijk[1]=j; ijk[2]=k;

        int orient[3];
        int ndim = -1;
        this->GetOrientationVector(rcvDataDescription,orient,ndim);

        int range[6];
        this->GetCellRefinedExtent(
              orient,ndim,i,j,k,nei.NeighborLevel,
              this->GetGridLevel(gridID),range);

        // Loop through the range of fine grid cells.
        int myIJK[3];
        for(int ii=IMIN(range); ii <= IMAX(range); ++ii)
660
        {
661
          for(int jj=JMIN(range); jj <= JMAX(range); ++jj)
662
          {
663
            for(int kk=KMIN(range); kk <= KMAX(range); ++kk)
664
            {
665 666
              myIJK[0]=ii; myIJK[1]=jj; myIJK[2]=kk;
              if( !this->IsNodeWithinExtent(ii,jj,kk,GhostedCellExtent))
667
              {
668
                continue;
669
              }
670 671 672

              if(this->IsNodeWithinExtent(i,j,k, NeighborCellExtent) &&
                 !this->IsNodeWithinExtent(ii,jj,kk,RegisteredGridCellExtent))
673
              {
674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692
                 // Sanity check!
                 assert("pre: RcvExtent is outside the GhostExtent!" &&
                    this->IsNodeWithinExtent(ii,jj,kk,GhostedGridExtent));
                 assert("pre: RcvExtent is outside the NeighborExtent" &&
                    this->IsNodeWithinExtent(i,j,k,NeighborExtent));

                 // Compute the source & target index
                 // Note: Since these indices are computed from a cell extent
                 // they correspond to a cell index.
                 vtkIdType sourceIdx =
                     vtkStructuredData::ComputePointIdForExtent(
                          NeighborCellExtent, ijk, this->DataDescription);

                 vtkIdType targetIdx =
                     vtkStructuredData::ComputePointIdForExtent(
                         GhostedCellExtent, myIJK, this->DataDescription);

                 if(this->CellCenteredDonorLevel[gridID][targetIdx] <
                     nei.NeighborLevel)
693
                 {
694 695 696 697 698
                   this->CopyFieldData(
                      this->GridCellData[nei.NeighborID],sourceIdx,
                      this->GhostedGridCellData[gridID],targetIdx);
                   this->CellCenteredDonorLevel[gridID][targetIdx]=
                       nei.NeighborLevel;
699 700
                 } // END if this is a finer solution
              } // END if
701

702 703 704
            } // END for all kk
          } // END for all jj
        } // END for all ii
705

706 707 708
      } // END for all k
    } // END for all j
  } // END for all i
709 710 711 712 713 714 715 716 717 718 719 720 721 722
}

//-----------------------------------------------------------------------------
void
vtkStructuredAMRGridConnectivity::GetLocalCellCentersFromFinerLevel(
    const int gridID, vtkStructuredAMRNeighbor &nei)
{
  assert( "pre: Expected a finer neighbor" &&
          (nei.NeighborLevel > nei.GridLevel) );

  // STEP 0: Get the grid's extent and cell extent
  int RegisteredGridExtent[6];
  this->GetGridExtent( gridID, RegisteredGridExtent );
  int RegisteredGridCellExtent[6];
723
  vtkStructuredData::GetCellExtentFromPointExtent(
724 725 726 727 728 729
      RegisteredGridExtent,RegisteredGridCellExtent,this->DataDescription);

  // STEP 1: Get the grid's ghosted extent and cell extent
  int GhostedGridExtent[6];
  this->GetGhostedExtent( gridID, GhostedGridExtent );
  int GhostedCellExtent[6];
730
  vtkStructuredData::GetCellExtentFromPointExtent(
731 732 733 734 735 736 737
      GhostedGridExtent,GhostedCellExtent,this->DataDescription);


  // STEP 3: Get the neighbor's extent and cell extent
  int NeighborExtent[6];
  this->GetGridExtent( nei.NeighborID, NeighborExtent );
  int NeighborCellExtent[6];
738
  vtkStructuredData::GetCellExtentFromPointExtent(
739 740 741
      NeighborExtent,NeighborCellExtent,this->DataDescription);

  // STEP 4: Get RcvCell extent
742 743 744
//  int rcvDataDescription =
//      vtkStructuredData::GetDataDescriptionFromExtent(
//                   const_cast<int*>(nei.RcvExtent));
745 746

  int RcvCellExtent[6];
747
  vtkStructuredData::GetCellExtentFromPointExtent(
748 749 750 751 752 753 754 755
      const_cast<int*>(nei.RcvExtent),RcvCellExtent);

  // STEP 5: Get receive node/cell extent w.r.t. this grid
  int GridRcvExtent[6];
 nei.GetReceiveExtentOnGrid(
     this->NumberOfGhostLayers,GhostedGridExtent,GridRcvExtent);
  int GridRcvDataDescription =
      vtkStructuredData::GetDataDescriptionFromExtent(GridRcvExtent);
756 757
//  assert("pre: mismatching data description" &&
//          (rcvDataDescription == GridRcvDataDescription) );
758
  int GridRcvCellExtent[6];
759
  vtkStructuredData::GetCellExtentFromPointExtent(
760 761 762 763
      GridRcvExtent,GridRcvCellExtent);

  int ijk[3];
  for(int i=IMIN(GridRcvCellExtent); i <= IMAX(GridRcvCellExtent); ++i)
764
  {
765
    for(int j=JMIN(GridRcvCellExtent); j <= JMAX(GridRcvCellExtent); ++j)
766
    {
767
      for(int k=KMIN(GridRcvCellExtent); k <= KMAX(GridRcvCellExtent); ++k)
768
      {
769 770 771
        ijk[0]=i; ijk[1]=j; ijk[2]=k;
        if(!this->IsNodeWithinExtent(i,j,k,RegisteredGridCellExtent) &&
            this->IsNodeWithinExtent(i,j,k,GhostedCellExtent))
772
        {
773 774 775 776 777 778 779 780
          // Compute target cell index. Note since a cell extent is given to
          // ComputePointIdForExtent, a cell index is returned.
          vtkIdType targetIdx =
              vtkStructuredData::ComputePointIdForExtent(
                  GhostedCellExtent,ijk,this->DataDescription);

          if(this->CellCenteredDonorLevel[gridID][targetIdx] <
             nei.NeighborLevel)
781
          {
782 783 784 785 786 787 788 789 790 791 792
            std::vector<vtkIdType> sourceIds;

            int range[6];
            int orient[3];
            int ndim;
            this->GetOrientationVector(GridRcvDataDescription,orient,ndim);
            this->GetCellRefinedExtent(
                orient,ndim,i,j,k,nei.GridLevel,nei.NeighborLevel,range);

            int rcvIJK[3];
            for(int ii=IMIN(range); ii <= IMAX(range); ++ii )
793
            {
794
              for(int jj=JMIN(range); jj <= JMAX(range); ++jj)
795
              {
796
                for(int kk=KMIN(range); kk <= KMAX(range); ++kk)
797
                {
798 799
                  rcvIJK[0]=ii; rcvIJK[1]=jj; rcvIJK[2]=kk;
                  if(this->IsNodeWithinExtent(ii,jj,kk,RcvCellExtent))
800
                  {
801 802 803 804
                    vtkIdType sourceIdx =
                        vtkStructuredData::ComputePointIdForExtent(
                            NeighborCellExtent,rcvIJK,this->DataDescription);
                    sourceIds.push_back(sourceIdx);
805 806 807 808
                  } // END if
                } // END for all kk
              } // END for all jj
            } // END for all ii
809

810
            if( !sourceIds.empty() )
811
            {
812 813
              this->AverageFieldData(
                  this->GridCellData[nei.NeighborID],&sourceIds[0],
814
                  static_cast<int>(sourceIds.size()),
815 816 817 818
                  this->GhostedGridCellData[gridID],targetIdx);

              this->CellCenteredDonorLevel[gridID][targetIdx]=
                  nei.NeighborLevel;
819
            }
820
            else
821
            {
822
              vtkWarningMacro("Empty list of sources!")
823
            }
824

825
          } // END if nei has finer resolution
826

827
        } // END if it is ghost cell
828

829 830 831
      } // END for all k
    } // END for all j
  } // END for all i
832 833 834 835 836 837 838 839 840 841 842 843
}

//-----------------------------------------------------------------------------
void
vtkStructuredAMRGridConnectivity::TransferLocalCellCenteredNeighborData(
    const int gridID, vtkStructuredAMRNeighbor &nei)
{
  int gridLevel = this->GetGridLevel( gridID );
  assert("pre: grid level mismatch!" && (gridLevel == nei.GridLevel) );

  // STEP 0: Check if the neighbor is strictly a child
  if( nei.RelationShip == vtkStructuredAMRNeighbor::CHILD )
844
  {
845 846 847
    // A child that is completely covered by this grid does not contribute to
    // its ghost-layers.
    return;
848
  }
849 850 851 852 853 854 855 856 857

  // STEP 1: Initialize the donor-level array if the array has not been
  // initialized before
  int GhostedGridExtent[6];
  this->GetGhostedExtent(gridID,GhostedGridExtent);
  int numCells = vtkStructuredData::GetNumberOfCells(
                       GhostedGridExtent,this->DataDescription);
  if( static_cast<int>(this->CellCenteredDonorLevel[gridID].size())
      != numCells)
858
  {
859
    this->CellCenteredDonorLevel[ gridID ].resize(numCells,-1);
860
  }
861 862 863

  // STEP 2: Fill data in the ghost levels
  if(gridLevel == nei.NeighborLevel)
864
  {
865
    this->GetLocalCellCentersAtSameLevel(gridID,nei);
866
  } // END if grids are at the same level
867
  else if(gridLevel < nei.NeighborLevel)
868
  {
869
    this->GetLocalCellCentersFromFinerLevel(gridID,nei);
870
  } // END if grid is coarser than the neighbor
871
  else
872
  {
873
    this->GetLocalCellCentersFromCoarserLevel(gridID,nei);
874
  } // END else grid is finer than the neighbor
875 876 877 878 879 880 881 882
}

//-----------------------------------------------------------------------------
void
vtkStructuredAMRGridConnectivity::TransferLocalNeighborData(
      const int gridID, vtkStructuredAMRNeighbor &nei)
{
  if( this->GetNodeCentered() )
883
  {
884
    this->TransferLocalNodeCenteredNeighborData(gridID,nei);
885
  }
886 887

  if( this->GetCellCentered() )
888
  {
889
    this->TransferLocalCellCenteredNeighborData(gridID,nei);
890
  }
891 892 893 894 895 896 897 898 899 900
}

//-----------------------------------------------------------------------------
void
vtkStructuredAMRGridConnectivity::TransferGhostDataFromNeighbors(
    const int gridID)
{
  // Sanity check
  assert("pre: gridID is out-of-bounds!" &&
       (gridID >= 0) && (gridID < static_cast<int>(this->NumberOfGrids)));
luz.paz's avatar
luz.paz committed
901
  assert("pre: Neighbors is not properly allocated" &&
902 903 904 905 906
       (this->NumberOfGrids==this->Neighbors.size() ) );

  this->CellCenteredDonorLevel.resize( this->NumberOfGrids );
  int NumNeis = static_cast<int>(this->Neighbors[gridID].size());
  for(int nei=0; nei < NumNeis; ++nei)
907
  {
908
    this->TransferLocalNeighborData(gridID,this->Neighbors[gridID][nei]);
909
  } // END for all neighbors
910 911 912 913 914 915 916 917
}

//-----------------------------------------------------------------------------
void
vtkStructuredAMRGridConnectivity::AverageFieldData(
    vtkFieldData *source, vtkIdType *sourceIds, const int N,
    vtkFieldData *target, vtkIdType targetIdx)
{
918 919 920
  assert( "pre: source field data is nullptr!" && (source != nullptr) );
  assert( "pre: target field data is nullptr!" && (target != nullptr) );
  assert( "pre: sourceIds is nullptr!" && (sourceIds != nullptr) );
921 922 923 924 925 926
  assert( "pre: N > 0" && (N > 0) );
  assert( "pre: source number of arrays does not match target!" &&
           source->GetNumberOfArrays()==target->GetNumberOfArrays() );

  int arrayIdx = 0;
  for( ; arrayIdx < source->GetNumberOfArrays(); ++arrayIdx )
927
  {
928 929
    // Get source array
    vtkDataArray *sourceArray = source->GetArray( arrayIdx );
930
    assert( "ERROR: encountered nullptr source array" && (sourceArray != nullptr) );
931 932 933

    // Get target array
    vtkDataArray *targetArray = target->GetArray( arrayIdx );
934
    assert( "ERROR: encountered nullptr target array" && (targetArray != nullptr) );
935 936 937 938 939 940 941 942 943 944 945 946 947 948 949

    // Sanity checks
    assert( "ERROR: target/source array name mismatch!" &&
      (strcmp( sourceArray->GetName(), targetArray->GetName() ) == 0 ) );
    assert( "ERROR: target/source array num components mismatch!" &&
      (sourceArray->GetNumberOfComponents()==
       targetArray->GetNumberOfComponents()));
    assert( "ERROR: targetIdx out-of-bounds!" &&
      (targetIdx >= 0) && (targetIdx < targetArray->GetNumberOfTuples() ) );

    int numComponents = sourceArray->GetNumberOfComponents();

    std::vector< double > averageTuple;
    averageTuple.resize(numComponents,0);
    for( int comp=0; comp < numComponents; ++comp)
950
    {
951
      for( int src=0; src < N; ++src )
952
      {
953 954 955 956 957
        vtkIdType sourceIdx = sourceIds[ src ];
        assert( "ERROR: sourceIdx out-of-bounds!" &&
         (sourceIdx >= 0) && (sourceIdx < sourceArray->GetNumberOfTuples()));

        averageTuple[comp] += sourceArray->GetComponent(sourceIdx,comp);
958
      }
959 960
      averageTuple[comp] /= N;
      targetArray->SetComponent(targetIdx,comp,averageTuple[comp]);
961
    } // END for all components
962

963
  } // END for all arrays
964 965 966 967 968 969 970 971
}

//-----------------------------------------------------------------------------
void
vtkStructuredAMRGridConnectivity::CopyFieldData(
    vtkFieldData *source, vtkIdType sourceIdx,
    vtkFieldData *target, vtkIdType targetIdx)
{
972 973
  assert( "pre: source field data is nullptr!" && (source != nullptr) );
  assert( "pre: target field data is nullptr!" && (target != nullptr) );
974 975 976 977 978
  assert( "pre: source number of arrays does not match target!" &&
          source->GetNumberOfArrays()==target->GetNumberOfArrays() );

  int arrayIdx = 0;
  for( ; arrayIdx < source->GetNumberOfArrays(); ++arrayIdx )
979
  {
980 981
    // Get source array
    vtkDataArray *sourceArray = source->GetArray( arrayIdx );
982
    assert( "ERROR: encountered nullptr source array" && (sourceArray != nullptr) );
983 984 985

    // Get target array
    vtkDataArray *targetArray = target->GetArray( arrayIdx );
986
    assert( "ERROR: encountered nullptr target array" && (targetArray != nullptr) );
987 988 989 990 991 992 993 994 995 996 997 998 999 1000

    // Sanity checks
    assert( "ERROR: target/source array name mismatch!" &&
      (strcmp( sourceArray->GetName(), targetArray->GetName() ) == 0 ) );
    assert( "ERROR: target/source array num components mismatch!" &&
      (sourceArray->GetNumberOfComponents()==
       targetArray->GetNumberOfComponents()));
    assert( "ERROR: sourceIdx out-of-bounds!" &&
      (sourceIdx >= 0) && (sourceIdx < sourceArray->GetNumberOfTuples() ) );
    assert( "ERROR: targetIdx out-of-bounds!" &&
      (targetIdx >= 0) && (targetIdx < targetArray->GetNumberOfTuples() ) );

    // Copy the tuple
    targetArray->SetTuple(targetIdx, sourceIdx, sourceArray);
1001
  } // END for all arrays
1002 1003 1004 1005 1006 1007 1008 1009 1010 1011
}

//-----------------------------------------------------------------------------
void
vtkStructuredAMRGridConnectivity::ComputeNeighborSendAndRcvExtent(
    const int gridID, const int N)
{
  // Sanity check
  assert( "pre: gridID is out-of-bounds!" &&
          (gridID >= 0) && (gridID < static_cast<int>(this->NumberOfGrids)));
luz.paz's avatar
luz.paz committed
1012
  assert( "pre: Neighbors is not properly allocated" &&
1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024
          (this->NumberOfGrids==this->Neighbors.size() ) );

  int gridRealExtent[6];
  this->GetGridExtent(gridID,gridRealExtent);

  int gridGhostedExtent[6];
  this->GetGhostedExtent(gridID,gridGhostedExtent);

  int neiRealExtent[6];

  int NumNeis = static_cast<int>(this->Neighbors[ gridID ].size());
  for(int nei=0; nei < NumNeis; ++nei)
1025
  {
1026 1027 1028 1029 1030
    this->GetGridExtent(
        this->Neighbors[gridID][nei].NeighborID,neiRealExtent);

    this->Neighbors[gridID][nei].ComputeSendAndReceiveExtent(
        gridRealExtent,gridGhostedExtent,neiRealExtent,this->WholeExtent,N);
1031
  } // END for all neighbors
1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045
}

//-----------------------------------------------------------------------------
void vtkStructuredAMRGridConnectivity::CreateGhostedMaskArrays(
        const int gridID)
{
  assert( "pre: gridID is out-of-bounds!" &&
       (gridID >= 0) && (gridID < static_cast<int>(this->NumberOfGrids)));
  assert( "pre: GhostedPointGhostArray has not been allocated" &&
       (this->NumberOfGrids == this->GhostedPointGhostArray.size()));
  assert( "pre: GhostedCellGhostArray has not been allocated" &&
       (this->NumberOfGrids == this->GhostedCellGhostArray.size()));

  // STEP 0: Initialize ghosted node and cell arrays
1046
  if( this->GhostedPointGhostArray[gridID] == nullptr )
1047
  {
1048
    this->GhostedPointGhostArray[ gridID ] = vtkUnsignedCharArray::New();
1049
  }
1050
  else
1051
  {
1052
    this->GhostedPointGhostArray[ gridID ]->Reset();
1053
  }
1054

1055
  if( this->GhostedCellGhostArray[ gridID ] == nullptr )
1056
  {
1057
    this->GhostedCellGhostArray[ gridID ] = vtkUnsignedCharArray::New();
1058
  }
1059
  else
1060
  {
1061
    this->GhostedCellGhostArray[ gridID ]->Reset();
1062
  }
1063 1064 1065 1066 1067 1068 1069

  // STEP 1: Get the ghosted extent
  int ghostExtent[6];
  this->GetGhostedExtent( gridID, ghostExtent );

  // STEP 2: Compute numNodes/numCells on the ghosted grid
  int numNodes =
1070
    vtkStructuredData::GetNumberOfPoints(ghostExtent,this->DataDescription);
1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088
  int numCells =
    vtkStructuredData::GetNumberOfCells(ghostExtent,this->DataDescription);

  // STEP 3: Allocate the ghosted node and cell arrays
  this->GhostedPointGhostArray[ gridID ]->Allocate( numNodes );
  this->GhostedCellGhostArray[ gridID ]->Allocate( numCells );

  // STEP 4: Get the registered extent of the grid
  int registeredGridExtent[6];
  this->GetGridExtent(gridID,registeredGridExtent);

  // STEP 5: Get normalized whole extent w.r.t. the level of this grid
  int normalizedWholeExt[6];
  this->GetWholeExtentAtLevel(this->GetGridLevel(gridID),normalizedWholeExt);

  // STEP 6: Fill ghosted points ghost array
  int ijk[3];
  for(int i=IMIN(ghostExtent); i <= IMAX(ghostExtent); ++i)
1089
  {
1090
    for(int j=JMIN(ghostExtent); j <= JMAX(ghostExtent); ++j)
1091
    {
1092
      for(int k=KMIN(ghostExtent); k <= KMAX(ghostExtent); ++k)
1093
      {
1094 1095 1096 1097 1098 1099
        ijk[0]=i; ijk[1]=j; ijk[2]=k;
        vtkIdType pntIdx =
           vtkStructuredData::ComputePointIdForExtent(
                   ghostExtent,ijk,this->DataDescription);

        if( this->IsNodeWithinExtent(i,j,k,registeredGridExtent))
1100
        {
1101 1102 1103 1104
          // Copy data from registered grid
          vtkIdType srcIdx =
             vtkStructuredData::ComputePointIdForExtent(
                 registeredGridExtent,ijk,this->DataDescription);
1105 1106
          unsigned char p = 0;
          if(this->GridPointGhostArrays[gridID])
1107
          {
1108
            p = this->GridPointGhostArrays[gridID]->GetValue(srcIdx);
1109
          }
1110
          this->GhostedPointGhostArray[ gridID ]->SetValue(pntIdx, p);
1111
        } // END if node within the registered extent
1112
        else
1113
        {
1114
          // The node is a ghost node
1115 1116
          unsigned char p = 0;
          p |= vtkDataSetAttributes::DUPLICATEPOINT;
1117
          if( this->IsNodeOnBoundaryOfExtent(i,j,k,normalizedWholeExt) )
1118
          {
1119 1120
            // We don't have BOUNDARY now but we might add
            // it in the future.
1121
          }
1122 1123

          this->GhostedPointGhostArray[ gridID ]->SetValue(pntIdx,p);
1124
        } // END else
1125

1126 1127 1128
      } // END for all k
    } // END for all j
  } // END for all i
1129 1130 1131

  // STEP 7: Fill ghosted cells ghost array
  int ghostCellExtent[6];
1132
  vtkStructuredData::GetCellExtentFromPointExtent(
1133 1134 1135
      ghostExtent,ghostCellExtent,this->DataDescription);

  int registeredCellExtent[6];
1136
  vtkStructuredData::GetCellExtentFromPointExtent(
1137 1138 1139
      registeredGridExtent,registeredCellExtent,this->DataDescription);

  for(int i=IMIN(ghostCellExtent); i <= IMAX(ghostCellExtent); ++i)
1140
  {
1141
    for(int j=JMIN(ghostCellExtent); j <= JMAX(ghostCellExtent); ++j)
1142
    {
1143
      for( int k=KMIN(ghostCellExtent); k <= KMAX(ghostCellExtent); ++k)
1144
      {
1145 1146 1147 1148 1149 1150
        ijk[0]=i; ijk[1]=j; ijk[2]=k;
        vtkIdType cellIdx =
            vtkStructuredData::ComputePointIdForExtent(
                ghostCellExtent,ijk,this->DataDescription);

        if( this->IsNodeWithinExtent(i,j,k,registeredCellExtent))
1151
        {
1152 1153 1154
          vtkIdType srcCellIdx =
              vtkStructuredData::ComputePointIdForExtent(
                  registeredCellExtent,ijk,this->DataDescription);
1155 1156
          unsigned char p = 0;
          if(this->GridCellGhostArrays[gridID])
1157
          {
1158
            p = this->GridCellGhostArrays[gridID]->GetValue(srcCellIdx);
1159
          }
1160 1161
          this->GhostedCellGhostArray[ gridID ]->SetValue(cellIdx, p);
        }
1162
        else
1163
        {
1164 1165
          unsigned char p = 0;
          p |= vtkDataSetAttributes::DUPLICATECELL;
1166
          this->GhostedCellGhostArray[ gridID ]->SetValue(cellIdx,p);
1167 1168 1169 1170
        }
      } // END for all k
    } // END for all j
  } // END for all i
1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187
}

//-----------------------------------------------------------------------------
void vtkStructuredAMRGridConnectivity::CreateGhostedExtent(
        const int gridId, const int N)
{
  assert("pre: gridId is out-of-bounds!" &&
         (gridId >= 0) && (gridId < static_cast<int>(this->NumberOfGrids)));
  assert("pre: ghosted extents vector has not been allocated!" &&
          (this->NumberOfGrids == this->GhostedExtents.size()/6 ) );
  assert("pre: number of ghost-layers requested should not be 0" &&
          (this->NumberOfGhostLayers > 0) );

  int ext[6];
  this->GetGridExtent( gridId, ext);

  switch( this->DataDescription )
1188
  {
1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200