Commit 91ada9f2 authored by George Zagaris's avatar George Zagaris
Browse files

WIP: Patch refinement

Change-Id: Iddce83442131f412446bba511cba641e3379797d
parent c60df9e8
......@@ -21,7 +21,9 @@
#include "vtkAMRUtilities.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkStructuredExtent.h"
#include "vtkCell.h"
#include "vtkMath.h"
#include <cassert>
......@@ -115,24 +117,84 @@ void vtkAMRGaussianPulseSource::GeneratePulseField(vtkUniformGrid* grid)
pulseField->Delete();
}
//------------------------------------------------------------------------------
vtkUniformGrid* vtkAMRGaussianPulseSource::RefinePatch(
vtkUniformGrid* parent, int patchExtent[6] )
{
assert("pre: parent grid is NULL!" && (parent!=NULL) );
int ext[6];
parent->GetExtent(ext);
assert("pre: patchExtent must be within the parent extent!" &&
vtkStructuredExtent::Smaller(patchExtent,ext));
double min[3];
double max[3];
double h[3];
double h0[3];
int ndim[3];
// Set some nominal values to ensure proper initialization
ndim[0] = ndim[1] = ndim[2] = 1;
min[0] = min[1] = min[2] =
max[0] = max[1] = max[2] = 0.0;
h[0] = h[1] = h[2] =
h0[0] = h0[1] = h0[2] = 0.5;
// STEP 0: Get min
int minIJK[3];
minIJK[0] = patchExtent[0];
minIJK[1] = patchExtent[2];
minIJK[2] = patchExtent[4];
vtkIdType minIdx = vtkStructuredData::ComputePointIdForExtent(ext,minIJK);
parent->GetPoint( minIdx, min );
// STEP 1: Get max
int maxIJK[3];
maxIJK[0] = patchExtent[1];
maxIJK[1] = patchExtent[3];
maxIJK[2] = patchExtent[5];
vtkIdType maxIdx = vtkStructuredData::ComputePointIdForExtent(ext,maxIJK);
parent->GetPoint( maxIdx, max );
// STEP 2: Compute the spacing of the refined patch and its dimensions
parent->GetSpacing( h0 );
for( int i=0; i < this->Dimension; ++i )
{
h[i] = h0[i]/static_cast<double>(this->RefinmentRatio);
ndim[i] = vtkMath::Floor(max[i]-min[i]/h[i]);
} // END for all dimensions
// STEP 3: Construct uniform grid for requested patch
vtkUniformGrid *grid = vtkUniformGrid::New();
grid->Initialize();
grid->SetOrigin( min );
grid->SetSpacing( h );
grid->SetDimensions(ndim);
// STEP 4: Compute Gaussian-Pulse field on patch
this->GeneratePulseField(grid);
return(grid);
}
//------------------------------------------------------------------------------
vtkUniformGrid* vtkAMRGaussianPulseSource::GetGrid(
double origin[3], double h[3], int ndim[3] )
{
vtkUniformGrid *grid = vtkUniformGrid::New();
grid->Initialize();
grid->SetOrigin( origin );
grid->SetSpacing( h );
grid->SetDimensions( ndim );
grid->SetOrigin(origin);
grid->SetSpacing(h);
grid->SetDimensions(ndim);
this->GeneratePulseField( grid );
return( grid );
this->GeneratePulseField(grid);
return(grid);
}
//------------------------------------------------------------------------------
void vtkAMRGaussianPulseSource::Generate2DDataSet( vtkOverlappingAMR *amr )
{
assert( "pre: input amr dataset is NULL" && (amr != NULL) );
assert( "pre: input amr dataset is NULL" && (amr!=NULL) );
int ndim[3];
double origin[3];
......@@ -140,6 +202,12 @@ void vtkAMRGaussianPulseSource::Generate2DDataSet( vtkOverlappingAMR *amr )
int blockId = 0;
int level = 0;
// Define the patches to be refined apriori
int patches[2][6] = {
{0,2,0,3,0,0},
{3,5,2,5,0,0}
};
// Root Block -- Block 0,0
ndim[0] = 6; ndim[1] = 6; ndim[2] = 1;
h[0] = h[1] = h[2] = this->RootSpacing[0];
......@@ -148,6 +216,17 @@ void vtkAMRGaussianPulseSource::Generate2DDataSet( vtkOverlappingAMR *amr )
level = 0;
vtkUniformGrid *grid = this->GetGrid(origin, h, ndim);
amr->SetDataSet(level,blockId,grid);
vtkUniformGrid *refinedPatch = NULL;
for( int patchIdx=0; patchIdx < 2; ++patchIdx )
{
refinedPatch = RefinePatch( grid, patches[patchIdx] );
assert("pre: refined grid is NULL" && (refinedPatch != NULL) );
amr->SetDataSet(level+1,patchIdx,refinedPatch);
refinedPatch->Delete();
refinedPatch = NULL;
}
grid->Delete();
vtkAMRUtilities::GenerateMetaData( amr, NULL );
amr->GenerateVisibilityArrays();
......@@ -173,6 +252,8 @@ void vtkAMRGaussianPulseSource::Generate3DDataSet( vtkOverlappingAMR *amr )
vtkUniformGrid *grid = this->GetGrid(origin, h, ndim);
amr->SetDataSet(level,blockId,grid);
grid->Delete();
vtkAMRUtilities::GenerateMetaData( amr, NULL );
amr->GenerateVisibilityArrays();
}
......
......@@ -98,9 +98,12 @@ class VTK_AMR_EXPORT vtkAMRGaussianPulseSource :
// Description:
// Computes the gaussian pulse at the given location based on the user
// supplied parameters for pulse width and origin.
double ComputePulseAt(const double x, const double y, const double z);
double ComputePulseAt( double pt[3] )
{return( this->ComputePulseAt(pt[0],pt[1],pt[2]) );}
double ComputePulseAt(const double x, const double y, const double z)
{
double xyz[3]; xyz[0]=x; xyz[1]=y; xyz[2]=z;
return( this->ComputePulseAt(xyz) );
}
double ComputePulseAt( double pt[3] );
// Description:
// Given the cell index w.r.t. to a uniform grid, this method computes the
......@@ -115,9 +118,13 @@ class VTK_AMR_EXPORT vtkAMRGaussianPulseSource :
// Description:
// Constructs a uniform grid path with the given origin/spacing and node
// dimensions.
// dimensions. The return grid serves as the root grid for the domain.
vtkUniformGrid* GetGrid( double origin[3], double h[3], int ndim[3] );
// Description:
// Constructs a refined patch from the given parent grid.
vtkUniformGrid* RefinePatch(vtkUniformGrid* parent, int patchExtent[6]);
// Description:
// Generate 2-D or 3-D DataSet
void Generate2DDataSet(vtkOverlappingAMR* amr);
......@@ -139,18 +146,18 @@ class VTK_AMR_EXPORT vtkAMRGaussianPulseSource :
//==============================================================================
// INLINE METHODS
//==============================================================================
inline double vtkAMRGaussianPulseSource::ComputePulseAt(
const double x, const double y, const double vtkNotUsed(z))
inline double vtkAMRGaussianPulseSource::ComputePulseAt(double pt[3])
{
double pulse = 0.0;
double r = 0.0;
double dx = x-this->PulseOrigin[0];
r += (dx*dx) / (this->PulseWidth[0]*this->PulseWidth[0]);
double dy = y-this->PulseOrigin[1];
r += (dy*dy) / (this->PulseWidth[1]*this->PulseWidth[1]);
for( int i=0; i < this->Dimension; ++i )
{
double d = pt[i]-this->PulseOrigin[i];
double d2 = d*d;
double L2 = this->PulseWidth[i]*this->PulseWidth[i];
r += d2/L2;
}
pulse = this->PulseAmplitude*std::exp( -r );
return( pulse );
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment