AMRDataTransferPipeline.cxx 8.49 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    AMRDataTransferPipeline.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.

=========================================================================*/
// .SECTION Description
//  A simple utility to demonstrate & test the parallel AMR functionality
//  and inter-block data transfer.

19
#include <cmath>
20
#include <cassert>
21
22
23
#include <mpi.h>

#include "vtkUniformGrid.h"
24
25
26
27
28
29
#include "vtkDataArray.h"
#include "vtkDoubleArray.h"
#include "vtkCell.h"
#include "vtkPoints.h"
#include "vtkPointData.h"
#include "vtkCellData.h"
30
#include "vtkAMRBox.h"
31
#include "vtkAMRConnectivityFilter.h"
George Zagaris's avatar
George Zagaris committed
32
#include "vtkAMRDataTransferFilter.h"
33
#include "vtkHierarchicalBoxDataSet.h"
34
#include "vtkXMLPHierarchicalBoxDataWriter.h"
35
#include "vtkMultiProcessController.h"
36
#include "vtkMPIController.h"
37
38
#include "vtkAMRUtilities.h"

39
40
41
42
43

// Global Variables
vtkMultiProcessController *Controller;

// Function Prototypes
44
45
vtkHierarchicalBoxDataSet* GetSerialAMRDataSet();
vtkHierarchicalBoxDataSet* GetParallelAMRDataSet();
46
vtkUniformGrid* GetGrid( double* origin,double* h,int* ndim );
47
void ComputeGaussianPulseError( vtkHierarchicalBoxDataSet *amrData );
48

49
50
51
//
// Main
//
52
53
int main( int argc, char **argv )
{
54
55
  Controller = vtkMPIController::New();
  Controller->Initialize( &argc, &argv );
56
  assert("pre: Controller != NULL" && (Controller != NULL ) );
57

58
  // STEP 0: Read In AMR data
59
  std::cout << "Constructing Sample AMR data!" << std::endl;
60
61
  std::cout.flush( );

62
63
64
65
66
67
68
69
70
71
72
73
  vtkHierarchicalBoxDataSet *amrData = NULL;

  if( Controller->GetNumberOfProcesses( ) == 3 )
    amrData = GetParallelAMRDataSet( );
  else if( Controller->GetNumberOfProcesses() == 1 )
    amrData = GetSerialAMRDataSet( );
  else
    {
      std::cerr << "Can only run with 1 or 3 MPI processes!\n";
      std::cerr.flush();
      return( -1 );
    }
74

75
76
  assert("pre: AMRData != NULL" && (amrData != NULL) );
  assert("pre: numLevels == 2" && (amrData->GetNumberOfLevels()==2) );
77
78
79
80
81
82
83
84
85

  std::cout << "Done reading!" << std::endl;
  std::cout.flush( );
  Controller->Barrier( );

  // STEP 1: Compute inter-block and inter-process connectivity
  std::cout << "Computing inter-block & inter-process connectivity!\n";
  std::cout.flush( );

86
  vtkAMRConnectivityFilter* connectivityFilter=vtkAMRConnectivityFilter::New( );
87
88
89
90
91
92
93
  connectivityFilter->SetController( Controller );
  connectivityFilter->SetAMRDataSet( amrData );
  connectivityFilter->ComputeConnectivity();

  std::cout << "Done computing connectivity!\n";
  std::cout.flush( );
  Controller->Barrier();
94

95
96
97
98
  connectivityFilter->Print( std::cout );
  Controller->Barrier();

  // STEP 2: Data transfer
George Zagaris's avatar
George Zagaris committed
99
100
101
102
103
104
105
  std::cout << "Transfering solution\n";
  std::cout.flush();

  vtkAMRDataTransferFilter* transferFilter = vtkAMRDataTransferFilter::New();

  transferFilter->SetController( Controller );
  transferFilter->SetAMRDataSet( amrData );
106
  transferFilter->SetNumberOfGhostLayers( 1 );
George Zagaris's avatar
George Zagaris committed
107
  transferFilter->SetRemoteConnectivity(
108
   connectivityFilter->GetRemoteConnectivity() );
George Zagaris's avatar
George Zagaris committed
109
  transferFilter->SetLocalConnectivity(
110
   connectivityFilter->GetLocalConnectivity() );
George Zagaris's avatar
George Zagaris committed
111
112
  transferFilter->Transfer();

113
114
115
116
  vtkHierarchicalBoxDataSet *newData = transferFilter->GetExtrudedData();
  assert( "extruded data is NULL!" && (newData != NULL) );
  ComputeGaussianPulseError( newData );

George Zagaris's avatar
George Zagaris committed
117
118
119
  std::cout << "[DONE]\n";
  std::cout.flush();
  Controller->Barrier();
120

121
122
123
124
125
  // STEP 3: CleanUp
  amrData->Delete();
  connectivityFilter->Delete();
  transferFilter->Delete();

126
  Controller->Finalize();
George Zagaris's avatar
George Zagaris committed
127
  Controller->Delete();
128
129
130
131
132
133
134
  return 0;
}

//=============================================================================
//                    Function Prototype Implementation
//=============================================================================

135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
//------------------------------------------------------------------------------
void ComputeGaussianPulseError( vtkHierarchicalBoxDataSet *data )
{
  // Sanity check
  assert( "pre: data is NULL!" && (data != NULL) );

  for( unsigned int level=0; level < data->GetNumberOfLevels(); ++level )
    {
      for( unsigned int idx=0; idx < data->GetNumberOfDataSets(level); ++idx)
        {

          vtkUniformGrid *ug = data->GetDataSet(level,idx);
          if( ug != NULL )
            {
              vtkDoubleArray *err = vtkDoubleArray::New();
              err->SetNumberOfComponents( 1 );
              err->SetNumberOfValues( ug->GetNumberOfCells() );

              for( int cell=0; cell < ug->GetNumberOfCells(); ++cell )
                {
                  //TODO: implement this
                }// END for all grid cells
            }

        } // END for all data
    } // END for all levels
}

163
//------------------------------------------------------------------------------
164
vtkHierarchicalBoxDataSet* GetSerialAMRDataSet( )
165
166
167
{

  vtkHierarchicalBoxDataSet *data = vtkHierarchicalBoxDataSet::New( );
168
169
170
171
172
173
174
175
  data->Initialize();

  double origin[3];
  double h[3];
  int    ndim[3];
  int    dim      =  2;
  int    blockId  = -1;
  int    level    = -1;
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
  int    rank     = 0;

  // BLOCK 0
  ndim[0]               = 25;
  ndim[1]               = 25;
  ndim[2]               = 1;
  origin[0]             = 0.0;
  origin[1]             = 0.0;
  origin[2]             = 0.0;
  blockId               = 0;
  level                 = 0;
  h[0] = h[1] = h[2]    = 0.5;
  vtkUniformGrid *grid0 = GetGrid( origin,h,ndim );
  data->SetDataSet( level,blockId,grid0);
  grid0->Delete();

  // BLOCK 1
  ndim[0]            = 11;
  ndim[1]            = 7;
  ndim[2]            = 1;
  origin[0]          = 1.5;
  origin[1]          = 1.5;
  origin[2]          = 0.0;
  blockId            = 0;
  level              = 1;
  h[0] = h[1] = h[2] = 0.25;
  vtkUniformGrid *grid1 = GetGrid( origin,h,ndim );
  data->SetDataSet( level,blockId,grid1);
  grid1->Delete();

  // BLOCK 2
  ndim[0]               = 11;
  ndim[1]               = 7;
  ndim[2]               = 1;
  origin[0]             = 1.0;
  origin[1]             = 3.0;
  origin[2]             = 0.0;
  blockId               = 1;
  level                 = 1;
  h[0] = h[1] = h[2]    = 0.25;
  vtkUniformGrid *grid2 = GetGrid( origin,h,ndim );
  data->SetDataSet( level,blockId,grid2);
  grid2->Delete();

  vtkAMRUtilities::GenerateMetaData( data, NULL );
  data->GenerateVisibilityArrays();
  return( data );
}

//------------------------------------------------------------------------------
vtkHierarchicalBoxDataSet* GetParallelAMRDataSet( )
{
  vtkHierarchicalBoxDataSet *data = vtkHierarchicalBoxDataSet::New( );
  data->Initialize();

  // TODO: imlpement this
232
233
234
235
236
237

  return( data );

}

//----------------------------------------------------------------------------
238
vtkUniformGrid* GetGrid( double *origin, double *spacing, int *ndim )
239
{
240
241
242
243
244
  vtkUniformGrid *grd = vtkUniformGrid::New();
  grd->Initialize();
  grd->SetOrigin( origin );
  grd->SetSpacing( spacing );
  grd->SetDimensions( ndim );
245

246
247
248
249
250
251
252
  // Apply Gaussian Pulse
  double x0 = 6.0;
  double y0 = 6.0;
  double lx = 12;
  double ly = 12;
  double a  = 0.1;

253
  vtkDoubleArray* xyz = vtkDoubleArray::New( );
254
  xyz->SetName( "GaussianPulse" );
255
256
257
258
259
  xyz->SetNumberOfComponents( 1 );
  xyz->SetNumberOfTuples( grd->GetNumberOfCells() );
  for( int cellIdx=0; cellIdx < grd->GetNumberOfCells(); ++cellIdx )
    {
      vtkCell* myCell = grd->GetCell( cellIdx);
260
      assert( "post: myCell != NULL" && (myCell != NULL) );
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277

      vtkPoints *cellPoints = myCell->GetPoints();

      double xyzCenter[3];
      xyzCenter[0] = 0.0;
      xyzCenter[1] = 0.0;
      xyzCenter[2] = 0.0;

      for( int cp=0; cp < cellPoints->GetNumberOfPoints(); ++cp )
        {
          double pnt[3];
          cellPoints->GetPoint( cp, pnt );
          xyzCenter[0] += pnt[0];
          xyzCenter[1] += pnt[1];
          xyzCenter[2] += pnt[2];
        } // END for all cell points

278
279
280
      double x = xyzCenter[0] / (cellPoints->GetNumberOfPoints());
      double y = xyzCenter[1] / (cellPoints->GetNumberOfPoints());
      double z = xyzCenter[2] / (cellPoints->GetNumberOfPoints());
281

282
283
284
      double r = ( ( (x-x0)*(x-x0) ) / ( lx*lx ) ) +
                 ( ( (y-y0)*(y-y0) ) / ( ly*ly ) );
      double f = a*std::exp( -r );
285
286
287
288
      xyz->SetTuple1(cellIdx,f);
    } // END for all cells

  grd->GetCellData()->AddArray(xyz);
George Zagaris's avatar
George Zagaris committed
289
  xyz->Delete();
290
  return grd;
291
}