diff --git a/IO/NetCDF/vtkNetCDFCAMReader.cxx b/IO/NetCDF/vtkNetCDFCAMReader.cxx
index 322076b239438cda0c93110cdf9249dc5d23f8cb..c6637d5a9ac087e6ae649ffe9ebb3a2d8f388bc4 100644
--- a/IO/NetCDF/vtkNetCDFCAMReader.cxx
+++ b/IO/NetCDF/vtkNetCDFCAMReader.cxx
@@ -18,6 +18,7 @@
#include "vtkDoubleArray.h"
#include "vtkFieldData.h"
#include "vtkFloatArray.h"
+#include "vtkIncrementalOctreePointLocator.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMath.h"
@@ -29,7 +30,7 @@
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkUnstructuredGrid.h"
-#include <map>
+#include <vector>
#include <vtk_netcdfcpp.h>
namespace
@@ -38,14 +39,30 @@ namespace
// a cell that wraps from the right side of the domain to the left side)
bool IsCellInverted(double points[4][3])
{
+ // We test the normal 3 points at a time. Not all grids are well-behaved
+ // i.e. consistenly use 0 or 360. We've had grid where 3 points on the left
+ // side, and just 1 on the right. Just checking the first 3 points (which is
+ // what ComputeNormal() does, we may (and do) miss a few cells.
+ // See BUG #0014897.
double normal[3];
- vtkPolygon::ComputeNormal(4, points[0], normal);
+ vtkPolygon::ComputeNormal(3, points[0], normal);
+ if(normal[2] > 0)
+ {
+ return true;
+ }
+ vtkPolygon::ComputeNormal(3, points[1], normal);
if(normal[2] > 0)
{
return true;
}
return false;
}
+
+ template <class T>
+ inline bool IsZero(const T& val)
+ {
+ return std::abs(val) < std::numeric_limits<T>::epsilon();
+ }
}
vtkStandardNewMacro(vtkNetCDFCAMReader);
@@ -60,7 +77,6 @@ vtkNetCDFCAMReader::vtkNetCDFCAMReader()
this->PointsFile = NULL;
this->ConnectivityFile = NULL;
this->SingleLevel = 0;
- this->CellLayerRight = 1;
this->TimeSteps = NULL;
this->NumberOfTimeSteps = 0;
this->SetNumberOfInputPorts(0);
@@ -330,6 +346,8 @@ int vtkNetCDFCAMReader::RequestData(
NcVar* lon = this->PointsFile->get_var("lon");
NcVar* lat = this->PointsFile->get_var("lat");
vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
+ output->SetPoints(points);
+
long numFilePoints = dimension->size();
if(lat == NULL || lon == NULL)
{
@@ -378,11 +396,19 @@ int vtkNetCDFCAMReader::RequestData(
// domain and only the points on the left boundary are included in
// the points file. if a cell uses a point that is on the left
// boundary and it should be on the right boundary we will have
- // to create that point. that's what boundaryPoints is used for.
- // note that if this->CellLayerRight is false then we do the opposite
- // and make the 'connecting' cell on the left side of the domain
- // we don't actually build the cells yet though.
- std::map<vtkIdType, vtkIdType> boundaryPoints;
+ // to create that point. That's what boundaryPoints is used for.
+ // The (index + numFilePoints) gives us the new point id, and the value
+ // for that in this array will correspond to the original point id that the
+ // boundaryPoint is duplicate of.
+ std::vector<vtkIdType> boundaryPoints;
+
+ // To avoid creating multiple duplicates, we we a
+ // vtkIncrementalOctreePointLocator.
+ vtkSmartPointer<vtkIncrementalOctreePointLocator> locator =
+ vtkSmartPointer<vtkIncrementalOctreePointLocator>::New();
+ locator->SetDataSet(output); // dataset only has points right now.
+ locator->BuildLocator();
+
dimension = this->ConnectivityFile->get_dim("ncells");
if(dimension == NULL)
{
@@ -411,8 +437,6 @@ int vtkNetCDFCAMReader::RequestData(
connectivity->set_cur(0, beginCell);
connectivity->get(&(cellConnectivity[0]), 4, numLocalCells);
- double leftSide = 1.;
- double rightSide = 359.;
for(long i=0;i<numLocalCells;i++)
{
vtkIdType pointIds[4];
@@ -422,45 +446,87 @@ int vtkNetCDFCAMReader::RequestData(
pointIds[j] = cellConnectivity[i+j*numLocalCells]-1;
points->GetPoint(pointIds[j], coords[j]);
}
- if(IsCellInverted(coords) == true)
+ if (IsCellInverted(coords) == true)
{
- for(int j=0;j<4;j++)
+ // First decide whether we're putting this cell on the 360 side (right) or on the
+ // 0 side (left). We decide this based on which side will have the
+ // smallest protrusion.
+ double delta = 0.0;
+ bool anchorLeft = false;
+ for (int j=0; j < 4; ++j)
+ {
+ // We're assured that coords[j][0] is in the range [0, 360].
+ // We just that fact to avoid having to do a std::abs() here.
+ double rightDelta = (360.0 - coords[j][0]);
+ double leftDelta = coords[j][0]; // i.e. (coords[j][0] - 0.0).
+ if (IsZero(rightDelta) || IsZero(leftDelta) || rightDelta == leftDelta)
+ {
+ // if the point is equidistant from both ends or is one of the ends,
+ // we let the other points in this cell dictate where the cell should
+ // anchor since this point can easily be anchored on either side with
+ // no side effects.
+ continue;
+ }
+ if (rightDelta < leftDelta)
+ {
+ if (rightDelta > delta)
+ {
+ delta = rightDelta;
+ anchorLeft = false;
+ }
+ }
+ else
+ {
+ if (leftDelta > delta)
+ {
+ delta = leftDelta;
+ anchorLeft = true;
+ }
+ }
+ }
+ // Once we've decided where we're anchoring we adjust the points.
+ for (int j=0; j < 4; ++j)
{
- if(this->CellLayerRight && coords[j][0] < leftSide)
+ if (anchorLeft)
{
- std::map<vtkIdType, vtkIdType>::iterator otherPoint =
- boundaryPoints.find(pointIds[j]);
- if(otherPoint != boundaryPoints.end())
- { // already made point on the right boundary
- cellConnectivity[i+j*numLocalCells] = otherPoint->second + 1;
+ // if coords[j] is closer to right (360), move it to the left.
+ if ( (360.0 - coords[j][0]) < coords[j][0] )
+ {
+ coords[j][0] -= 360.0;
}
else
- { // need to make point on the right boundary
- vtkIdType index =
- points->InsertNextPoint(coords[j][0]+360., coords[j][1], coords[j][2]);
- cellConnectivity[i+j*numLocalCells] = index+1;
- boundaryPoints[pointIds[j]] = index;
+ {
+ continue;
}
}
- else if(this->CellLayerRight == 0 && coords[j][0] > rightSide)
+ else
{
- std::map<vtkIdType, vtkIdType>::iterator otherPoint =
- boundaryPoints.find(pointIds[j]);
- if(otherPoint != boundaryPoints.end())
- { // already made point on the right boundary
- cellConnectivity[i+j*numLocalCells] = otherPoint->second+1;
+ // if coords[j] is closer to left (0), move it to the right
+ if ( coords[j][0] < (360.0 - coords[j][0]) )
+ {
+ coords[j][0] += 360.0;
}
else
- { // need to make point on the right boundary
- vtkIdType index =
- points->InsertNextPoint(coords[j][0]-360., coords[j][1], coords[j][2]);
- cellConnectivity[i+j*numLocalCells] = index+1;
- boundaryPoints[pointIds[j]] = index;
+ {
+ continue;
}
}
+ // Okay, we have moved the coords. Update the boundaryPoints so which
+ // original point id is this new point id a clone of.
+ vtkIdType newPtId;
+ if (locator->InsertUniquePoint(coords[j], newPtId) == 1)
+ {
+ // if a new point was indeed inserted, we need to update the
+ // boundaryPoints to keep track of it.
+ assert(newPtId >= numFilePoints && pointIds[j] < newPtId);
+ assert(static_cast<vtkIdType>(boundaryPoints.size()) == (newPtId-numFilePoints));
+ boundaryPoints.push_back(pointIds[j]);
+ }
+ cellConnectivity[i+j*numLocalCells] = (newPtId + 1); // note: 1-indexed.
}
}
}
+ locator = NULL; // release the locator memory.
// we now have all of the points at a single level. build them up
// for the rest of the levels before creating the cells.
@@ -484,7 +550,7 @@ int vtkNetCDFCAMReader::RequestData(
}
points->Modified();
- output->SetPoints(points);
+ points->Squeeze();
this->SetProgress(.5); // educated guess for progress
@@ -601,13 +667,14 @@ int vtkNetCDFCAMReader::RequestData(
pointData->CopyAllOn();
pointData->CopyAllocate(output->GetPointData(),
output->GetNumberOfPoints());
- for(std::map<vtkIdType, vtkIdType>::const_iterator it=
- boundaryPoints.begin();it!=boundaryPoints.end();it++)
+
+ vtkIdType newPtId = numFilePoints;
+ for (std::vector<vtkIdType>::const_iterator it=
+ boundaryPoints.begin(); it!=boundaryPoints.end(); ++it, ++newPtId)
{
for(long lev=0;lev<numLocalCellLevels+1-this->SingleLevel;lev++)
{
- pointData->CopyData(pointData, it->first+lev*numPointsPerLevel,
- it->second+lev*numPointsPerLevel);
+ pointData->CopyData(pointData, (*it), newPtId);
}
}
@@ -755,6 +822,23 @@ bool vtkNetCDFCAMReader::GetPartitioning(
return true;
}
+//----------------------------------------------------------------------------
+#if !defined(VTK_LEGACY_REMOVE)
+void vtkNetCDFCAMReader::SetCellLayerRight(int)
+{
+ VTK_LEGACY_BODY(vtkNetCDFCAMReader::SetCellLayerRight, "VTK 6.3");
+}
+#endif
+
+//----------------------------------------------------------------------------
+#if !defined(VTK_LEGACY_REMOVE)
+int vtkNetCDFCAMReader::GetCellLayerRight()
+{
+ VTK_LEGACY_BODY(vtkNetCDFCAMReader::GetCellLayerRight, "VTK 6.3");
+ return 0;
+}
+#endif
+
//----------------------------------------------------------------------------
void vtkNetCDFCAMReader::PrintSelf(ostream& os, vtkIndent indent)
{
@@ -765,7 +849,6 @@ void vtkNetCDFCAMReader::PrintSelf(ostream& os, vtkIndent indent)
(this->ConnectivityFileName ? this->ConnectivityFileName : "(NULL)")
<< endl;
os << indent << "SingleLevel: " << this->SingleLevel << endl;
- os << indent << "CellLayerRight: " << this->CellLayerRight << endl;
if(this->PointsFile)
{
os << indent << "PointsFile: " << this->PointsFile << endl;
diff --git a/IO/NetCDF/vtkNetCDFCAMReader.h b/IO/NetCDF/vtkNetCDFCAMReader.h
index 6100dfc0c015208fa91203902f357b6020fc5208..43f0b637e53e2cd97d702cf3d817a9791c2e1e89 100644
--- a/IO/NetCDF/vtkNetCDFCAMReader.h
+++ b/IO/NetCDF/vtkNetCDFCAMReader.h
@@ -67,8 +67,10 @@ public:
// Specify which "side" of the domain to add the connecting
// cells at. 0 indicates left side and 1 indicates right side.
// The default is the right side.
- vtkSetMacro(CellLayerRight, int);
- vtkGetMacro(CellLayerRight, int);
+ // @deprecated This method is no longer supported. The reader automatically
+ // decides which side to pad cells on. Using this method has no effect.
+ VTK_LEGACY(void SetCellLayerRight(int));
+ VTK_LEGACY(int GetCellLayerRight());
protected:
vtkNetCDFCAMReader();
@@ -110,8 +112,6 @@ private:
int SingleLevel;
- int CellLayerRight;
-
double * TimeSteps;
long NumberOfTimeSteps;