Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-16 17:22:55 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Add specialization of topology map fetch for regular point coords

Browse Source 
In the special case where you are loading the point coordinates for a
structured grid in a point to cell map (an important use case), create a
VecRectilinearPointCoordinates rather than build a Vec of the values.
This will activate the cell specalizations in previous commits.

These changes also added some flat-to-logical index conversion and vice
versa in ConnectivityStructuredInternals. This change also fixed a bug
in getting cells attached to points in 2D grids. (Actually, technically
someone else fixed it and checked it in first. The changes were merged
during a rebase.)

I also added a specalization to Vec for 1D that implicitly converts
between the 1D Vec and the component. This can be convenient when
templating on the Vec length.
This commit is contained in:
Kenneth Moreland 2015-09-02 11:08:25 -07:00
parent b58543297a
commit 08f9c04fab
6 changed files with 526 additions and 122 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
vtkm/Types.h
View File

@ -980,6 +980,27 @@ public:
}
};
// Vectors of size 1 should implicitly convert between the scalar and the
// vector. Otherwise, it should behave the same.
template<typename T>
class Vec<T,1> : public detail::VecBase<T, 1, Vec<T,1> >
{
typedef detail::VecBase<T, 1, Vec<T,1> > Superclass;
public:
VTKM_EXEC_CONT_EXPORT Vec() {}
VTKM_EXEC_CONT_EXPORT Vec(const T& value) : Superclass(value) { }
template<typename OtherType>
VTKM_EXEC_CONT_EXPORT Vec(const Vec<OtherType, 1> &src) : Superclass(src) { }
VTKM_EXEC_CONT_EXPORT
operator T() const
{
return this->Components[0];
}
};
//-----------------------------------------------------------------------------
// Specializations for common tuple sizes (with special names).

28
vtkm/exec/ConnectivityStructured.h
View File

@ -85,6 +85,34 @@ public:
return Helper::GetIndices(this->Internals,index);
}
VTKM_EXEC_CONT_EXPORT
vtkm::Vec<vtkm::Id,Dimension>
FlatToLogicalPointIndex(vtkm::Id flatPointIndex) const
{
return this->Internals.FlatToLogicalPointIndex(flatPointIndex);
}
VTKM_EXEC_CONT_EXPORT
vtkm::Id LogicalToFlatPointIndex(
const vtkm::Vec<vtkm::Id,Dimension> &logicalPointIndex) const
{
return this->Internals.LogicalToFlatPointIndex(logicalPointIndex);
}
VTKM_EXEC_CONT_EXPORT
vtkm::Vec<vtkm::Id,Dimension>
FlatToLogicalCellIndex(vtkm::Id flatCellIndex) const
{
return this->Internals.FlatToLogicalCellIndex(flatCellIndex);
}
VTKM_EXEC_CONT_EXPORT
vtkm::Id LogicalToFlatCellIndex(
const vtkm::Vec<vtkm::Id,Dimension> &logicalCellIndex) const
{
return this->Internals.LogicalToFlatCellIndex(logicalCellIndex);
}
private:
InternalsType Internals;
};

78
vtkm/exec/arg/FetchTagArrayTopologyMapIn.h
View File

@ -23,6 +23,13 @@
#include <vtkm/exec/arg/AspectTagDefault.h>
#include <vtkm/exec/arg/Fetch.h>
#include <vtkm/TopologyElementTag.h>
#include <vtkm/internal/ArrayPortalUniformPointCoordinates.h>
#include <vtkm/exec/ConnectivityStructured.h>
#include <vtkm/VecRectilinearPointCoordinates.h>
#include <vtkm/exec/internal/VecFromPortalPermute.h>
namespace vtkm {
@ -69,6 +76,77 @@ struct FetchArrayTopologyMapInImplementation
}
};
VTKM_EXEC_EXPORT
vtkm::VecRectilinearPointCoordinates<1>
make_VecRectilinearPointCoordinates(
const vtkm::Vec<vtkm::FloatDefault,3> &origin,
const vtkm::Vec<vtkm::FloatDefault,3> &spacing,
const vtkm::Vec<vtkm::Id,1> &logicalId)
{
vtkm::Vec<vtkm::FloatDefault,3> offsetOrigin(
origin[0] + spacing[0]*static_cast<vtkm::FloatDefault>(logicalId[0]),
origin[1],
origin[2]);
return vtkm::VecRectilinearPointCoordinates<1>(offsetOrigin, spacing);
}
VTKM_EXEC_EXPORT
vtkm::VecRectilinearPointCoordinates<2>
make_VecRectilinearPointCoordinates(
const vtkm::Vec<vtkm::FloatDefault,3> &origin,
const vtkm::Vec<vtkm::FloatDefault,3> &spacing,
const vtkm::Vec<vtkm::Id,2> &logicalId)
{
vtkm::Vec<vtkm::FloatDefault,3> offsetOrigin(
origin[0] + spacing[0]*static_cast<vtkm::FloatDefault>(logicalId[0]),
origin[1] + spacing[1]*static_cast<vtkm::FloatDefault>(logicalId[1]),
origin[2]);
return vtkm::VecRectilinearPointCoordinates<2>(offsetOrigin, spacing);
}
VTKM_EXEC_EXPORT
vtkm::VecRectilinearPointCoordinates<3>
make_VecRectilinearPointCoordinates(
const vtkm::Vec<vtkm::FloatDefault,3> &origin,
const vtkm::Vec<vtkm::FloatDefault,3> &spacing,
const vtkm::Vec<vtkm::Id,3> &logicalId)
{
vtkm::Vec<vtkm::FloatDefault,3> offsetOrigin(
origin[0] + spacing[0]*static_cast<vtkm::FloatDefault>(logicalId[0]),
origin[1] + spacing[1]*static_cast<vtkm::FloatDefault>(logicalId[1]),
origin[2] + spacing[2]*static_cast<vtkm::FloatDefault>(logicalId[2]));
return vtkm::VecRectilinearPointCoordinates<3>(offsetOrigin, spacing);
}
template<vtkm::IdComponent NumDimensions>
struct FetchArrayTopologyMapInImplementation<
vtkm::exec::ConnectivityStructured<vtkm::TopologyElementTagPoint,
vtkm::TopologyElementTagCell,
NumDimensions>,
vtkm::internal::ArrayPortalUniformPointCoordinates>
{
typedef vtkm::exec::ConnectivityStructured<vtkm::TopologyElementTagPoint,
vtkm::TopologyElementTagCell,
NumDimensions> ConnectivityType;
typedef vtkm::VecRectilinearPointCoordinates<NumDimensions> ValueType;
VTKM_EXEC_EXPORT
static ValueType Load(
vtkm::Id index,
const ConnectivityType &connectivity,
const vtkm::internal::ArrayPortalUniformPointCoordinates &field)
{
// This works because the logical cell index is the same as the logical
// point index of the first point on the cell.
return vtkm::exec::arg::detail::make_VecRectilinearPointCoordinates(
field.GetOrigin(),
field.GetSpacing(),
connectivity.FlatToLogicalCellIndex(index));
}
};
} // namespace detail
template<typename Invocation, vtkm::IdComponent ParameterIndex>

1
vtkm/exec/arg/testing/CMakeLists.txt
View File

@ -23,6 +23,7 @@ set(unit_tests
UnitTestFetchArrayDirectIn.cxx
UnitTestFetchArrayDirectInOut.cxx
UnitTestFetchArrayDirectOut.cxx
UnitTestFetchArrayTopologyMapIn.cxx
UnitTestFetchExecObject.cxx
UnitTestFetchWorkIndex.cxx
)

233
vtkm/exec/arg/testing/UnitTestFetchArrayTopologyMapIn.cxx Normal file
View File

@ -0,0 +1,233 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt 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.
//
// Copyright 2014 Sandia Corporation.
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#include <vtkm/exec/arg/FetchTagArrayTopologyMapIn.h>
#include <vtkm/internal/FunctionInterface.h>
#include <vtkm/internal/Invocation.h>
#include <vtkm/testing/Testing.h>
namespace {
static const vtkm::Id ARRAY_SIZE = 10;
template<typename T>
struct TestPortal
{
typedef T ValueType;
VTKM_EXEC_CONT_EXPORT
vtkm::Id GetNumberOfValues() const { return ARRAY_SIZE; }
VTKM_EXEC_CONT_EXPORT
ValueType Get(vtkm::Id index) const {
VTKM_TEST_ASSERT(index >= 0, "Bad portal index.");
VTKM_TEST_ASSERT(index < this->GetNumberOfValues(), "Bad portal index.");
return TestValue(index, ValueType());
}
};
struct NullParam { };
template<vtkm::IdComponent InputDomainIndex,
vtkm::IdComponent ParamIndex,
typename T>
struct FetchArrayTopologyMapInTests
{
template<typename Invocation>
void TryInvocation(const Invocation &invocation) const
{
typedef vtkm::exec::arg::Fetch<
vtkm::exec::arg::FetchTagArrayTopologyMapIn,
vtkm::exec::arg::AspectTagDefault,
Invocation,
ParamIndex> FetchType;
FetchType fetch;
typename FetchType::ValueType value = fetch.Load(0, invocation);
VTKM_TEST_ASSERT(value.GetNumberOfComponents() == 8,
"Topology fetch got wrong number of components.");
VTKM_TEST_ASSERT(test_equal(value[0], TestValue(0, T())),
"Got invalid value from Load.");
VTKM_TEST_ASSERT(test_equal(value[1], TestValue(1, T())),
"Got invalid value from Load.");
VTKM_TEST_ASSERT(test_equal(value[2], TestValue(3, T())),
"Got invalid value from Load.");
VTKM_TEST_ASSERT(test_equal(value[3], TestValue(2, T())),
"Got invalid value from Load.");
VTKM_TEST_ASSERT(test_equal(value[4], TestValue(4, T())),
"Got invalid value from Load.");
VTKM_TEST_ASSERT(test_equal(value[5], TestValue(5, T())),
"Got invalid value from Load.");
VTKM_TEST_ASSERT(test_equal(value[6], TestValue(7, T())),
"Got invalid value from Load.");
VTKM_TEST_ASSERT(test_equal(value[7], TestValue(6, T())),
"Got invalid value from Load.");
}
void operator()() const
{
std::cout << "Trying ArrayTopologyMapIn fetch on parameter " << ParamIndex
<< " with type " << vtkm::testing::TypeName<T>::Name()
<< std::endl;
typedef vtkm::internal::FunctionInterface<
void(NullParam,NullParam,NullParam,NullParam,NullParam)>
BaseFunctionInterface;
vtkm::internal::ConnectivityStructuredInternals<3> connectivityInternals;
connectivityInternals.SetPointDimensions(vtkm::Id3(2,2,2));
vtkm::exec::ConnectivityStructured<
vtkm::TopologyElementTagPoint,vtkm::TopologyElementTagCell,3>
connectivity(connectivityInternals);
this->TryInvocation(vtkm::internal::make_Invocation<InputDomainIndex>(
BaseFunctionInterface()
.Replace<InputDomainIndex>(connectivity)
.template Replace<ParamIndex>(TestPortal<T>()),
NullParam(),
NullParam()));
}
};
struct TryType
{
template<typename T>
void operator()(T) const
{
FetchArrayTopologyMapInTests<3,1,T>()();
FetchArrayTopologyMapInTests<1,2,T>()();
FetchArrayTopologyMapInTests<2,3,T>()();
FetchArrayTopologyMapInTests<1,4,T>()();
FetchArrayTopologyMapInTests<1,5,T>()();
}
};
template<vtkm::IdComponent NumDimensions,
vtkm::IdComponent ParamIndex,
typename Invocation>
void TryStructuredPointCoordinatesInvocation(const Invocation &invocation)
{
vtkm::exec::arg::Fetch<
vtkm::exec::arg::FetchTagArrayTopologyMapIn,
vtkm::exec::arg::AspectTagDefault,
Invocation,
ParamIndex> fetch;
vtkm::Vec<vtkm::FloatDefault,3> origin =
TestValue(0, vtkm::Vec<vtkm::FloatDefault,3>());
vtkm::Vec<vtkm::FloatDefault,3> spacing =
TestValue(1, vtkm::Vec<vtkm::FloatDefault,3>());
vtkm::VecRectilinearPointCoordinates<NumDimensions> value =
fetch.Load(0, invocation);
VTKM_TEST_ASSERT(test_equal(value.GetOrigin(), origin), "Bad origin.");
VTKM_TEST_ASSERT(test_equal(value.GetSpacing(), spacing), "Bad spacing.");
origin[0] += spacing[0];
value = fetch.Load(1, invocation);
VTKM_TEST_ASSERT(test_equal(value.GetOrigin(), origin), "Bad origin.");
VTKM_TEST_ASSERT(test_equal(value.GetSpacing(), spacing), "Bad spacing.");
}
template<vtkm::IdComponent NumDimensions>
void TryStructuredPointCoordinates(
const vtkm::exec::ConnectivityStructured<
vtkm::TopologyElementTagPoint,vtkm::TopologyElementTagCell,NumDimensions> &connectivity,
const vtkm::internal::ArrayPortalUniformPointCoordinates &coordinates)
{
typedef vtkm::internal::FunctionInterface<
void(NullParam,NullParam,NullParam,NullParam,NullParam)>
BaseFunctionInterface;
// Try with topology in argument 1 and point coordinates in argument 2
TryStructuredPointCoordinatesInvocation<NumDimensions,2>(
vtkm::internal::make_Invocation<1>(
BaseFunctionInterface()
.Replace<1>(connectivity)
.template Replace<2>(coordinates),
NullParam(),
NullParam())
);
// Try again with topology in argument 3 and point coordinates in argument 1
TryStructuredPointCoordinatesInvocation<NumDimensions,1>(
vtkm::internal::make_Invocation<3>(
BaseFunctionInterface()
.Replace<3>(connectivity)
.template Replace<1>(coordinates),
NullParam(),
NullParam())
);
}
void TryStructuredPointCoordinates()
{
std::cout << "*** Fetching special case of uniform point coordinates. *****"
<< std::endl;
vtkm::internal::ArrayPortalUniformPointCoordinates coordinates(
vtkm::Id3(3,2,2),
TestValue(0, vtkm::Vec<vtkm::FloatDefault,3>()),
TestValue(1, vtkm::Vec<vtkm::FloatDefault,3>()));
std::cout << "3D" << std::endl;
vtkm::internal::ConnectivityStructuredInternals<3> connectivityInternals3d;
connectivityInternals3d.SetPointDimensions(vtkm::Id3(3,2,2));
vtkm::exec::ConnectivityStructured<
vtkm::TopologyElementTagPoint, vtkm::TopologyElementTagCell, 3>
connectivity3d(connectivityInternals3d);
TryStructuredPointCoordinates(connectivity3d, coordinates);
std::cout << "2D" << std::endl;
vtkm::internal::ConnectivityStructuredInternals<2> connectivityInternals2d;
connectivityInternals2d.SetPointDimensions(vtkm::Id2(3,2));
vtkm::exec::ConnectivityStructured<
vtkm::TopologyElementTagPoint, vtkm::TopologyElementTagCell, 2>
connectivity2d(connectivityInternals2d);
TryStructuredPointCoordinates(connectivity2d, coordinates);
std::cout << "1D" << std::endl;
vtkm::internal::ConnectivityStructuredInternals<1> connectivityInternals1d;
connectivityInternals1d.SetPointDimensions(3);
vtkm::exec::ConnectivityStructured<
vtkm::TopologyElementTagPoint, vtkm::TopologyElementTagCell, 1>
connectivity1d(connectivityInternals1d);
TryStructuredPointCoordinates(connectivity1d, coordinates);
}
void TestArrayTopologyMapIn()
{
vtkm::testing::Testing::TryTypes(TryType(),
vtkm::TypeListTagCommon());
TryStructuredPointCoordinates();
}
} // anonymous namespace
int UnitTestFetchArrayTopologyMapIn(int, char *[])
{
return vtkm::testing::Testing::Run(TestArrayTopologyMapIn);
}

287
vtkm/internal/ConnectivityStructuredInternals.h
View File

@ -121,6 +121,30 @@ public:
return cellIds;
}
VTKM_EXEC_CONT_EXPORT
vtkm::Id FlatToLogicalPointIndex(vtkm::Id flatPointIndex) const
{
return flatPointIndex;
}
VTKM_EXEC_CONT_EXPORT
vtkm::Id LogicalToFlatPointIndex(vtkm::Id logicalPointIndex) const
{
return logicalPointIndex;
}
VTKM_EXEC_CONT_EXPORT
vtkm::Id FlatToLogicalCellIndex(vtkm::Id flatCellIndex) const
{
return flatCellIndex;
}
VTKM_EXEC_CONT_EXPORT
vtkm::Id LogicalToFlatCellIndex(vtkm::Id logicalCellIndex) const
{
return logicalCellIndex;
}
VTKM_CONT_EXPORT
void PrintSummary(std::ostream &out) const
{
@ -187,13 +211,13 @@ public:
typedef vtkm::CellShapeTagQuad CellShapeTag;
VTKM_EXEC_CONT_EXPORT
vtkm::Vec<vtkm::Id,NUM_POINTS_IN_CELL> GetPointsOfCell(vtkm::Id index) const
vtkm::Vec<vtkm::Id,NUM_POINTS_IN_CELL>
GetPointsOfCell(vtkm::Id cellIndex) const
{
vtkm::Id i, j;
this->CalculateLogicalCellIndices(index, i, j);
vtkm::Id2 ij = this->FlatToLogicalCellIndex(cellIndex);
vtkm::Vec<vtkm::Id,NUM_POINTS_IN_CELL> pointIds;
pointIds[0] = j*this->PointDimensions[0] + i;
pointIds[0] = this->LogicalToFlatPointIndex(ij);
pointIds[1] = pointIds[0] + 1;
pointIds[2] = pointIds[1] + this->PointDimensions[0];
pointIds[3] = pointIds[2] - 1;
@ -203,44 +227,76 @@ public:
VTKM_EXEC_CONT_EXPORT
vtkm::IdComponent GetNumberOfCellsIncidentOnPoint(vtkm::Id pointIndex) const
{
vtkm::Id i, j;
this->CalculateLogicalPointIndices(pointIndex, i, j);
vtkm::Id2 ij = this->FlatToLogicalPointIndex(pointIndex);
return
(static_cast<vtkm::IdComponent>((i > 0) && (j > 0))
+ static_cast<vtkm::IdComponent>((i < this->PointDimensions[0]-1) && (j > 0))
+ static_cast<vtkm::IdComponent>((i > 0) && (j < this->PointDimensions[1]-1))
(static_cast<vtkm::IdComponent>((ij[0] > 0) && (ij[1] > 0))
+ static_cast<vtkm::IdComponent>((ij[0] < this->PointDimensions[0]-1) && (ij[1] > 0))
+ static_cast<vtkm::IdComponent>((ij[0] > 0) && (ij[1] < this->PointDimensions[1]-1))
+ static_cast<vtkm::IdComponent>(
(i < this->PointDimensions[0]-1) && (j < this->PointDimensions[1]-1)));
(ij[0] < this->PointDimensions[0]-1) && (ij[1] < this->PointDimensions[1]-1)));
}
VTKM_EXEC_CONT_EXPORT
vtkm::VecVariable<vtkm::Id,MAX_CELL_TO_POINT>
GetCellsOfPoint(vtkm::Id index) const
GetCellsOfPoint(vtkm::Id pointIndex) const
{
vtkm::VecVariable<vtkm::Id,MAX_CELL_TO_POINT> cellIds;
vtkm::Id i, j;
this->CalculateLogicalPointIndices(index, i, j);
if ((i > 0) && (j > 0))
vtkm::Id2 ij = this->FlatToLogicalPointIndex(pointIndex);
if ((ij[0] > 0) && (ij[1] > 0))
{
cellIds.Append(this->CalculateCellIndex(i-1, j-1));
cellIds.Append(this->LogicalToFlatCellIndex(ij - vtkm::Id2(1,1)));
}
if ((i < this->PointDimensions[0]-1) && (j > 0))
if ((ij[0] < this->PointDimensions[0]-1) && (ij[1] > 0))
{
cellIds.Append(this->CalculateCellIndex(i , j-1));
cellIds.Append(this->LogicalToFlatCellIndex(ij - vtkm::Id2(0,1)));
}
if ((i > 0) && (j < this->PointDimensions[1]-1))
if ((ij[0] > 0) && (ij[1] < this->PointDimensions[1]-1))
{
cellIds.Append(this->CalculateCellIndex(i-1, j ));
cellIds.Append(this->LogicalToFlatCellIndex(ij - vtkm::Id2(1,0)));
}
if ((i < this->PointDimensions[0]-1) && (j < this->PointDimensions[1]-1))
if ((ij[0] < this->PointDimensions[0]-1) && (ij[1] < this->PointDimensions[1]-1))
{
cellIds.Append(this->CalculateCellIndex(i , j ));
cellIds.Append(this->LogicalToFlatCellIndex(ij));
}
return cellIds;
}
VTKM_EXEC_CONT_EXPORT
vtkm::Id2 FlatToLogicalPointIndex(vtkm::Id flatPointIndex) const
{
vtkm::Id2 logicalPointIndex;
logicalPointIndex[0] = flatPointIndex % this->PointDimensions[0];
logicalPointIndex[1] = flatPointIndex / this->PointDimensions[0];
return logicalPointIndex;
}
VTKM_EXEC_CONT_EXPORT
vtkm::Id LogicalToFlatPointIndex(const vtkm::Id2 &logicalPointIndex) const
{
return logicalPointIndex[0]
+ this->PointDimensions[0]*logicalPointIndex[1];
}
VTKM_EXEC_CONT_EXPORT
vtkm::Id2 FlatToLogicalCellIndex(vtkm::Id flatCellIndex) const
{
vtkm::Id2 cellDimensions = this->GetCellDimensions();
vtkm::Id2 logicalCellIndex;
logicalCellIndex[0] = flatCellIndex % cellDimensions[0];
logicalCellIndex[1] = flatCellIndex / cellDimensions[0];
return logicalCellIndex;
}
VTKM_EXEC_CONT_EXPORT
vtkm::Id LogicalToFlatCellIndex(const vtkm::Id2 &logicalCellIndex) const
{
vtkm::Id2 cellDimensions = this->GetCellDimensions();
return logicalCellIndex[0]
+ cellDimensions[0]*logicalCellIndex[1];
}
VTKM_CONT_EXPORT
void PrintSummary(std::ostream &out) const
{
@ -251,28 +307,6 @@ public:
private:
vtkm::Id2 PointDimensions;
private:
VTKM_EXEC_CONT_EXPORT
void CalculateLogicalPointIndices(vtkm::Id index, vtkm::Id &i, vtkm::Id &j) const
{
vtkm::Id2 pointDimensions = this->GetPointDimensions();
i = index % pointDimensions[0];
j = index / pointDimensions[0];
}
VTKM_EXEC_CONT_EXPORT
void CalculateLogicalCellIndices(vtkm::Id index, vtkm::Id &i, vtkm::Id &j) const
{
vtkm::Id2 cellDimensions = this->GetCellDimensions();
i = index % cellDimensions[0];
j = index / cellDimensions[0];
}
VTKM_EXEC_CONT_EXPORT
vtkm::Id CalculateCellIndex(vtkm::Id i, vtkm::Id j) const
{
vtkm::Id2 cellDimensions = this->GetCellDimensions();
return j*cellDimensions[0] + i;
}
};
//3 D specialization.
@ -332,13 +366,12 @@ public:
typedef vtkm::CellShapeTagHexahedron CellShapeTag;
VTKM_EXEC_CONT_EXPORT
vtkm::Vec<vtkm::Id,NUM_POINTS_IN_CELL> GetPointsOfCell(vtkm::Id index) const
vtkm::Vec<vtkm::Id,NUM_POINTS_IN_CELL> GetPointsOfCell(vtkm::Id cellIndex) const
{
vtkm::Id i, j, k;
this->CalculateLogicalCellIndices(index, i, j, k);
vtkm::Id3 ijk = this->FlatToLogicalCellIndex(cellIndex);
vtkm::Vec<vtkm::Id,NUM_POINTS_IN_CELL> pointIds;
pointIds[0] = (k * this->PointDimensions[1] + j) * this->PointDimensions[0] + i;
pointIds[0] = (ijk[2] * this->PointDimensions[1] + ijk[1]) * this->PointDimensions[0] + ijk[0];
pointIds[1] = pointIds[0] + 1;
pointIds[2] = pointIds[1] + this->PointDimensions[0];
pointIds[3] = pointIds[2] - 1;
@ -353,79 +386,76 @@ public:
VTKM_EXEC_CONT_EXPORT
vtkm::IdComponent GetNumberOfCellsIncidentOnPoint(vtkm::Id pointIndex) const
{
vtkm::Id i, j, k;
this->CalculateLogicalPointIndices(pointIndex, i, j, k);
vtkm::Id3 ijk = this->FlatToLogicalPointIndex(pointIndex);
return (
static_cast<vtkm::IdComponent>((i > 0) && (j > 0) && (k > 0))
+ static_cast<vtkm::IdComponent>((i < this->PointDimensions[0]-1) && (j > 0) && (k > 0))
+ static_cast<vtkm::IdComponent>((i > 0) && (j < this->PointDimensions[1]-1) && (k > 0))
+ static_cast<vtkm::IdComponent>((i < this->PointDimensions[0]-1) &&
(j < this->PointDimensions[1]-1) &&
(k > 0))
+ static_cast<vtkm::IdComponent>((i > 0) && (j > 0) && (k < this->PointDimensions[2]-1))
+ static_cast<vtkm::IdComponent>((i < this->PointDimensions[0]-1) &&
(j > 0) &&
(k < this->PointDimensions[2]-1))
+ static_cast<vtkm::IdComponent>((i > 0) &&
(j < this->PointDimensions[1]-1) &&
(k < this->PointDimensions[2]-1))
+ static_cast<vtkm::IdComponent>((i < this->PointDimensions[0]-1) &&
(j < this->PointDimensions[1]-1) &&
(k < this->PointDimensions[2]-1))
static_cast<vtkm::IdComponent>((ijk[0] > 0) && (ijk[1] > 0) && (ijk[2] > 0))
+ static_cast<vtkm::IdComponent>((ijk[0] < this->PointDimensions[0]-1) && (ijk[1] > 0) && (ijk[2] > 0))
+ static_cast<vtkm::IdComponent>((ijk[0] > 0) && (ijk[1] < this->PointDimensions[1]-1) && (ijk[2] > 0))
+ static_cast<vtkm::IdComponent>((ijk[0] < this->PointDimensions[0]-1) &&
(ijk[1] < this->PointDimensions[1]-1) &&
(ijk[2] > 0))
+ static_cast<vtkm::IdComponent>((ijk[0] > 0) && (ijk[1] > 0) && (ijk[2] < this->PointDimensions[2]-1))
+ static_cast<vtkm::IdComponent>((ijk[0] < this->PointDimensions[0]-1) &&
(ijk[1] > 0) &&
(ijk[2] < this->PointDimensions[2]-1))
+ static_cast<vtkm::IdComponent>((ijk[0] > 0) &&
(ijk[1] < this->PointDimensions[1]-1) &&
(ijk[2] < this->PointDimensions[2]-1))
+ static_cast<vtkm::IdComponent>((ijk[0] < this->PointDimensions[0]-1) &&
(ijk[1] < this->PointDimensions[1]-1) &&
(ijk[2] < this->PointDimensions[2]-1))
);
}
VTKM_EXEC_CONT_EXPORT
vtkm::VecVariable<vtkm::Id,MAX_CELL_TO_POINT>
GetCellsOfPoint(vtkm::Id index) const
GetCellsOfPoint(vtkm::Id pointIndex) const
{
vtkm::VecVariable<vtkm::Id,MAX_CELL_TO_POINT> cellIds;
vtkm::Id i, j, k;
vtkm::Id3 ijk = this->FlatToLogicalPointIndex(pointIndex);
this->CalculateLogicalPointIndices(index, i, j, k);
if ((i > 0) && (j > 0) && (k > 0))
if ((ijk[0] > 0) && (ijk[1] > 0) && (ijk[2] > 0))
{
cellIds.Append(this->CalculateCellIndex(i-1, j-1, k-1));
cellIds.Append(this->LogicalToFlatCellIndex(ijk - vtkm::Id3(1,1,1)));
}
if ((i < this->PointDimensions[0]-1) && (j > 0) && (k > 0))
if ((ijk[0] < this->PointDimensions[0]-1) && (ijk[1] > 0) && (ijk[2] > 0))
{
cellIds.Append(this->CalculateCellIndex(i , j-1, k-1));
cellIds.Append(this->LogicalToFlatCellIndex(ijk - vtkm::Id3(0,1,1)));
}
if ((i > 0) && (j < this->PointDimensions[1]-1) && (k > 0))
if ((ijk[0] > 0) && (ijk[1] < this->PointDimensions[1]-1) && (ijk[2] > 0))
{
cellIds.Append(this->CalculateCellIndex(i-1, j , k-1));
cellIds.Append(this->LogicalToFlatCellIndex(ijk - vtkm::Id3(1,0,1)));
}
if ((i < this->PointDimensions[0]-1) &&
(j < this->PointDimensions[1]-1) &&
(k > 0))
if ((ijk[0] < this->PointDimensions[0]-1) &&
(ijk[1] < this->PointDimensions[1]-1) &&
(ijk[2] > 0))
{
cellIds.Append(this->CalculateCellIndex(i , j , k-1));
cellIds.Append(this->LogicalToFlatCellIndex(ijk - vtkm::Id3(0,0,1)));
}
if ((i > 0) && (j > 0) && (k < this->PointDimensions[2]-1))
if ((ijk[0] > 0) && (ijk[1] > 0) && (ijk[2] < this->PointDimensions[2]-1))
{
cellIds.Append(this->CalculateCellIndex(i-1, j-1, k));
cellIds.Append(this->LogicalToFlatCellIndex(ijk - vtkm::Id3(1,1,0)));
}
if ((i < this->PointDimensions[0]-1) &&
(j > 0) &&
(k < this->PointDimensions[2]-1))
if ((ijk[0] < this->PointDimensions[0]-1) &&
(ijk[1] > 0) &&
(ijk[2] < this->PointDimensions[2]-1))
{
cellIds.Append(this->CalculateCellIndex(i , j-1, k));
cellIds.Append(this->LogicalToFlatCellIndex(ijk - vtkm::Id3(0,1,0)));
}
if ((i > 0) &&
(j < this->PointDimensions[1]-1) &&
(k < this->PointDimensions[2]-1))
if ((ijk[0] > 0) &&
(ijk[1] < this->PointDimensions[1]-1) &&
(ijk[2] < this->PointDimensions[2]-1))
{
cellIds.Append(this->CalculateCellIndex(i-1, j , k));
cellIds.Append(this->LogicalToFlatCellIndex(ijk - vtkm::Id3(1,0,0)));
}
if ((i < this->PointDimensions[0]-1) &&
(j < this->PointDimensions[1]-1) &&
(k < this->PointDimensions[2]-1))
if ((ijk[0] < this->PointDimensions[0]-1) &&
(ijk[1] < this->PointDimensions[1]-1) &&
(ijk[2] < this->PointDimensions[2]-1))
{
cellIds.Append(this->CalculateCellIndex(i , j , k));
cellIds.Append(this->LogicalToFlatCellIndex(ijk));
}
return cellIds;
@ -439,38 +469,51 @@ public:
out<<std::endl;
}
VTKM_EXEC_CONT_EXPORT
vtkm::Id3 FlatToLogicalPointIndex(vtkm::Id flatPointIndex) const
{
vtkm::Id3 logicalPointIndex;
vtkm::Id pointDims01 = this->PointDimensions[0] * this->PointDimensions[1];
logicalPointIndex[2] = flatPointIndex / pointDims01;
vtkm::Id indexij = flatPointIndex % pointDims01;
logicalPointIndex[1] = indexij / this->PointDimensions[0];
logicalPointIndex[0] = indexij % this->PointDimensions[0];
return logicalPointIndex;
}
VTKM_EXEC_CONT_EXPORT
vtkm::Id LogicalToFlatPointIndex(const vtkm::Id3 &logicalPointIndex) const
{
return logicalPointIndex[0]
+ this->PointDimensions[0]*(logicalPointIndex[1]
+ this->PointDimensions[1]*logicalPointIndex[2]);
}
VTKM_EXEC_CONT_EXPORT
vtkm::Id3 FlatToLogicalCellIndex(vtkm::Id flatCellIndex) const
{
vtkm::Id3 cellDimensions = this->GetCellDimensions();
vtkm::Id3 logicalCellIndex;
vtkm::Id cellDims01 = cellDimensions[0] * cellDimensions[1];
logicalCellIndex[2] = flatCellIndex / cellDims01;
vtkm::Id indexij = flatCellIndex % cellDims01;
logicalCellIndex[1] = indexij / cellDimensions[0];
logicalCellIndex[0] = indexij % cellDimensions[0];
return logicalCellIndex;
}
VTKM_EXEC_CONT_EXPORT
vtkm::Id LogicalToFlatCellIndex(const vtkm::Id3 &logicalCellIndex) const
{
vtkm::Id3 cellDimensions = this->GetCellDimensions();
return logicalCellIndex[0]
+ cellDimensions[0]*(logicalCellIndex[1]
+ cellDimensions[1]*logicalCellIndex[2]);
}
private:
vtkm::Id3 PointDimensions;
private:
VTKM_EXEC_CONT_EXPORT
void CalculateLogicalCellIndices(vtkm::Id index, vtkm::Id &i, vtkm::Id &j, vtkm::Id &k) const
{
vtkm::Id3 cellDimensions = this->GetCellDimensions();
vtkm::Id cellDims01 = cellDimensions[0] * cellDimensions[1];
k = index / cellDims01;
vtkm::Id indexij = index % cellDims01;
j = indexij / cellDimensions[0];
i = indexij % cellDimensions[0];
}
VTKM_EXEC_CONT_EXPORT
void CalculateLogicalPointIndices(vtkm::Id index, vtkm::Id &i, vtkm::Id &j, vtkm::Id &k) const
{
vtkm::Id pointDims01 = this->PointDimensions[0] * this->PointDimensions[1];
k = index / pointDims01;
vtkm::Id indexij = index % pointDims01;
j = indexij / this->PointDimensions[0];
i = indexij % this->PointDimensions[0];
}
VTKM_EXEC_CONT_EXPORT
vtkm::Id CalculateCellIndex(vtkm::Id i, vtkm::Id j, vtkm::Id k) const
{
vtkm::Id3 cellDimensions = this->GetCellDimensions();
return (k * cellDimensions[1] + j) * cellDimensions[0] + i;
}
};
// We may want to generalize this class depending on how ConnectivityExplicit