AMRDataTransferPipeline.cxx 10.3 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 <sstream>
21
#include <cassert>
22
23
24
#include <mpi.h>

#include "vtkUniformGrid.h"
25
26
27
28
29
30
#include "vtkDataArray.h"
#include "vtkDoubleArray.h"
#include "vtkCell.h"
#include "vtkPoints.h"
#include "vtkPointData.h"
#include "vtkCellData.h"
31
#include "vtkAMRBox.h"
32
#include "vtkAMRConnectivityFilter.h"
George Zagaris's avatar
George Zagaris committed
33
#include "vtkAMRDataTransferFilter.h"
34
#include "vtkHierarchicalBoxDataSet.h"
35
#include "vtkXMLHierarchicalBoxDataWriter.h"
36
#include "vtkMultiProcessController.h"
37
#include "vtkMPIController.h"
38
#include "vtkAMRUtilities.h"
39
#include "vtkAMRDualMeshExtractor.h"
40
41
42
43
44

// Global Variables
vtkMultiProcessController *Controller;

// Function Prototypes
45
void WriteAMRData( vtkHierarchicalBoxDataSet *amrData, std::string prefix );
46
47
vtkHierarchicalBoxDataSet* GetSerialAMRDataSet();
vtkHierarchicalBoxDataSet* GetParallelAMRDataSet();
48
vtkUniformGrid* GetGrid( double* origin,double* h,int* ndim );
49
void ComputeGaussianPulseError( vtkHierarchicalBoxDataSet *amrData );
50
double ComputePulseAt( vtkUniformGrid *grid, const int cellIdx );
51

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

61
  // STEP 0: Read In AMR data
62
  std::cout << "Constructing Sample AMR data!" << std::endl;
63
64
  std::cout.flush( );

65
66
67
68
69
70
71
72
73
74
75
76
  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 );
    }
77

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

  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( );

89
  vtkAMRConnectivityFilter* connectivityFilter=vtkAMRConnectivityFilter::New( );
90
91
92
93
94
95
96
  connectivityFilter->SetController( Controller );
  connectivityFilter->SetAMRDataSet( amrData );
  connectivityFilter->ComputeConnectivity();

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

98
99
100
101
  connectivityFilter->Print( std::cout );
  Controller->Barrier();

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

  vtkAMRDataTransferFilter* transferFilter = vtkAMRDataTransferFilter::New();

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

116
117
118
  vtkHierarchicalBoxDataSet *newData = transferFilter->GetExtrudedData();
  assert( "extruded data is NULL!" && (newData != NULL) );
  ComputeGaussianPulseError( newData );
119
  WriteAMRData( newData, "NEWDATA" );
120

George Zagaris's avatar
George Zagaris committed
121
122
123
  std::cout << "[DONE]\n";
  std::cout.flush();
  Controller->Barrier();
124

125
126
127
128
129
  vtkAMRDualMeshExtractor *dme = vtkAMRDualMeshExtractor::New();
  dme->SetInput( newData );
  dme->Update();
  dme->WriteMultiBlockData( dme->GetOutput() );

130
  // STEP 3: CleanUp
131
  dme->Delete();
132
133
134
135
  amrData->Delete();
  connectivityFilter->Delete();
  transferFilter->Delete();

136
  Controller->Finalize();
George Zagaris's avatar
George Zagaris committed
137
  Controller->Delete();
138
139
140
141
142
143
144
  return 0;
}

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

145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
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
196
197
198
199
200
201
202
203
204
205
//------------------------------------------------------------------------------
void WriteAMRData( vtkHierarchicalBoxDataSet *amrData, std::string prefix )
{
  assert( "pre: AMR Data is NULL!" && (amrData != NULL) );

  vtkXMLHierarchicalBoxDataWriter *myAMRWriter=
    vtkXMLHierarchicalBoxDataWriter::New();

  std::ostringstream oss;
  oss << prefix << "." << myAMRWriter->GetDefaultFileExtension();

  myAMRWriter->SetFileName( oss.str().c_str() );
  myAMRWriter->SetInput( amrData );
  myAMRWriter->Write();
  myAMRWriter->Delete();
}

//------------------------------------------------------------------------------
double ComputePulseAt( vtkUniformGrid *grid, const int cellIdx )
{
  // Sanity check
  assert( "pre: grid != NULL" && (grid != NULL) );
  assert( "pre: cellIdx in bounds" &&
         ( (cellIdx >= 0) && (cellIdx < grid->GetNumberOfCells()) ) );

  // Gaussian Pulse parameters
  double x0 = 6.0;
  double y0 = 6.0;
  double lx = 12;
  double ly = 12;
  double a  = 0.1;

  vtkCell* myCell = grid->GetCell( cellIdx);
  assert( "post: myCell != NULL" && (myCell != NULL) );

  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

  double x = xyzCenter[0] / (cellPoints->GetNumberOfPoints());
  double y = xyzCenter[1] / (cellPoints->GetNumberOfPoints());
  double z = xyzCenter[2] / (cellPoints->GetNumberOfPoints());
  double r = ( ( (x-x0)*(x-x0) ) / ( lx*lx ) ) +
             ( ( (y-y0)*(y-y0) ) / ( ly*ly ) );
  double f = a*std::exp( -r );

  return( f );
}

206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
//------------------------------------------------------------------------------
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();
221
              err->SetName( "err" );
222
223
224
              err->SetNumberOfComponents( 1 );
              err->SetNumberOfValues( ug->GetNumberOfCells() );

225
226
227
228
229
230
231
              vtkCellData *CD = ug->GetCellData();
              assert( "pre: Grid must have cell data!" &&
                      (CD->GetNumberOfArrays() > 0) );

              vtkDataArray *pulse = CD->GetArray( "GaussianPulse" );
              assert( "pre: Pulse data not found on grid!" && (pulse!=NULL) );

232
233
              for( int cell=0; cell < ug->GetNumberOfCells(); ++cell )
                {
234
235
236
237
                  double expected = ComputePulseAt( ug, cell );
                  double actual   = pulse->GetComponent( cell, 0 );
                  double abserr   = std::abs( (actual-expected) );
                  err->SetComponent( cell, 0, abserr );
238
239
240
241
                } // END for all grid cells

              CD->AddArray( err );
              err->Delete();
242
243
244
245
246
247
            }

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

248
//------------------------------------------------------------------------------
249
vtkHierarchicalBoxDataSet* GetSerialAMRDataSet( )
250
251
252
{

  vtkHierarchicalBoxDataSet *data = vtkHierarchicalBoxDataSet::New( );
253
254
255
256
257
258
259
260
  data->Initialize();

  double origin[3];
  double h[3];
  int    ndim[3];
  int    dim      =  2;
  int    blockId  = -1;
  int    level    = -1;
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
  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
317
318
319
320
321
322

  return( data );

}

//----------------------------------------------------------------------------
323
vtkUniformGrid* GetGrid( double *origin, double *spacing, int *ndim )
324
{
325
326
327
328
329
  vtkUniformGrid *grd = vtkUniformGrid::New();
  grd->Initialize();
  grd->SetOrigin( origin );
  grd->SetSpacing( spacing );
  grd->SetDimensions( ndim );
330
331

  vtkDoubleArray* xyz = vtkDoubleArray::New( );
332
  xyz->SetName( "GaussianPulse" );
333
334
335
336
  xyz->SetNumberOfComponents( 1 );
  xyz->SetNumberOfTuples( grd->GetNumberOfCells() );
  for( int cellIdx=0; cellIdx < grd->GetNumberOfCells(); ++cellIdx )
    {
337
      xyz->SetTuple1(cellIdx, ComputePulseAt(grd,cellIdx) );
338
339
340
    } // END for all cells

  grd->GetCellData()->AddArray(xyz);
George Zagaris's avatar
George Zagaris committed
341
  xyz->Delete();
342
  return grd;
343
}