33edcc940d
There are several switch statements in ParametricCoordinates.h that have a default clause that is an error condition. The code was simply calling the RaiseError method on the worklet before quiting. Although the behavior is supposed to be similar to throwing an exception, flow actually does still continue. These error default clauses kept going without breaking from the switch or initializing the variables the method is supposed to set. Although the values should not matter once the error is raised, I was getting compiler warnings about it. It's also bad practice as certain NaN values can cause the program to signal out, so it is best to set the data to something.
1021 lines
36 KiB
C++
1021 lines
36 KiB
C++
//============================================================================
|
|
// 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 2015 Sandia Corporation.
|
|
// Copyright 2015 UT-Battelle, LLC.
|
|
// Copyright 2015 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.
|
|
//============================================================================
|
|
#ifndef vtk_m_exec_ParametricCoordinates_h
|
|
#define vtk_m_exec_ParametricCoordinates_h
|
|
|
|
#include <vtkm/CellShape.h>
|
|
#include <vtkm/Math.h>
|
|
#include <vtkm/exec/Assert.h>
|
|
#include <vtkm/exec/CellDerivative.h>
|
|
#include <vtkm/exec/CellInterpolate.h>
|
|
#include <vtkm/exec/FunctorBase.h>
|
|
#include <vtkm/exec/NewtonsMethod.h>
|
|
|
|
namespace vtkm {
|
|
namespace exec {
|
|
|
|
//-----------------------------------------------------------------------------
|
|
template<typename ParametricCoordType>
|
|
VTKM_EXEC_EXPORT
|
|
void ParametricCoordinatesCenter(vtkm::IdComponent numPoints,
|
|
vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagEmpty,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
VTKM_ASSERT_EXEC(numPoints == 0, worklet);
|
|
pcoords[0] = 0;
|
|
pcoords[1] = 0;
|
|
pcoords[2] = 0;
|
|
}
|
|
|
|
template<typename ParametricCoordType>
|
|
VTKM_EXEC_EXPORT
|
|
void ParametricCoordinatesCenter(vtkm::IdComponent numPoints,
|
|
vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagVertex,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
VTKM_ASSERT_EXEC(numPoints == 1, worklet);
|
|
pcoords[0] = 0;
|
|
pcoords[1] = 0;
|
|
pcoords[2] = 0;
|
|
}
|
|
|
|
template<typename ParametricCoordType>
|
|
VTKM_EXEC_EXPORT
|
|
void ParametricCoordinatesCenter(vtkm::IdComponent numPoints,
|
|
vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagLine,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
VTKM_ASSERT_EXEC(numPoints == 2, worklet);
|
|
pcoords[0] = 0.5;
|
|
pcoords[1] = 0;
|
|
pcoords[2] = 0;
|
|
}
|
|
|
|
template<typename ParametricCoordType>
|
|
VTKM_EXEC_EXPORT
|
|
void ParametricCoordinatesCenter(vtkm::IdComponent numPoints,
|
|
vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagTriangle,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
VTKM_ASSERT_EXEC(numPoints == 3, worklet);
|
|
pcoords[0] = static_cast<ParametricCoordType>(1.0/3.0);
|
|
pcoords[1] = static_cast<ParametricCoordType>(1.0/3.0);
|
|
pcoords[2] = 0;
|
|
}
|
|
|
|
template<typename ParametricCoordType>
|
|
VTKM_EXEC_EXPORT
|
|
void ParametricCoordinatesCenter(vtkm::IdComponent numPoints,
|
|
vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagPolygon,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
VTKM_ASSERT_EXEC(numPoints > 0, worklet);
|
|
switch (numPoints)
|
|
{
|
|
case 1:
|
|
ParametricCoordinatesCenter(
|
|
numPoints, pcoords, vtkm::CellShapeTagVertex(), worklet);
|
|
break;
|
|
case 2:
|
|
ParametricCoordinatesCenter(
|
|
numPoints, pcoords, vtkm::CellShapeTagLine(), worklet);
|
|
break;
|
|
case 3:
|
|
ParametricCoordinatesCenter(
|
|
numPoints, pcoords, vtkm::CellShapeTagTriangle(), worklet);
|
|
break;
|
|
default:
|
|
pcoords[0] = 0.5;
|
|
pcoords[1] = 0.5;
|
|
pcoords[2] = 0;
|
|
break;
|
|
}
|
|
}
|
|
|
|
template<typename ParametricCoordType>
|
|
VTKM_EXEC_EXPORT
|
|
void ParametricCoordinatesCenter(vtkm::IdComponent numPoints,
|
|
vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagQuad,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
VTKM_ASSERT_EXEC(numPoints == 4, worklet);
|
|
pcoords[0] = 0.5;
|
|
pcoords[1] = 0.5;
|
|
pcoords[2] = 0;
|
|
}
|
|
|
|
template<typename ParametricCoordType>
|
|
VTKM_EXEC_EXPORT
|
|
void ParametricCoordinatesCenter(vtkm::IdComponent numPoints,
|
|
vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagTetra,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
VTKM_ASSERT_EXEC(numPoints == 4, worklet);
|
|
pcoords[0] = 0.25;
|
|
pcoords[1] = 0.25;
|
|
pcoords[2] = 0.25;
|
|
}
|
|
|
|
template<typename ParametricCoordType>
|
|
VTKM_EXEC_EXPORT
|
|
void ParametricCoordinatesCenter(vtkm::IdComponent numPoints,
|
|
vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagHexahedron,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
VTKM_ASSERT_EXEC(numPoints == 8, worklet);
|
|
pcoords[0] = 0.5;
|
|
pcoords[1] = 0.5;
|
|
pcoords[2] = 0.5;
|
|
}
|
|
|
|
template<typename ParametricCoordType>
|
|
VTKM_EXEC_EXPORT
|
|
void ParametricCoordinatesCenter(vtkm::IdComponent numPoints,
|
|
vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagWedge,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
VTKM_ASSERT_EXEC(numPoints == 6, worklet);
|
|
pcoords[0] = static_cast<ParametricCoordType>(1.0/3.0);
|
|
pcoords[1] = static_cast<ParametricCoordType>(1.0/3.0);
|
|
pcoords[2] = 0.5;
|
|
}
|
|
|
|
template<typename ParametricCoordType>
|
|
VTKM_EXEC_EXPORT
|
|
void ParametricCoordinatesCenter(vtkm::IdComponent numPoints,
|
|
vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagPyramid,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
VTKM_ASSERT_EXEC(numPoints == 5, worklet);
|
|
pcoords[0] = 0.5;
|
|
pcoords[1] = 0.5;
|
|
pcoords[2] = static_cast<ParametricCoordType>(0.2);
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
/// Returns the parametric center of the given cell shape with the given number
|
|
/// of points.
|
|
///
|
|
template<typename ParametricCoordType>
|
|
VTKM_EXEC_EXPORT
|
|
void ParametricCoordinatesCenter(vtkm::IdComponent numPoints,
|
|
vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagGeneric shape,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
switch (shape.Id)
|
|
{
|
|
vtkmGenericCellShapeMacro(ParametricCoordinatesCenter(numPoints,
|
|
pcoords,
|
|
CellShapeTag(),
|
|
worklet));
|
|
default:
|
|
worklet.RaiseError("Bad shape given to ParametricCoordinatesCenter.");
|
|
pcoords[0] = pcoords[1] = pcoords[2] = 0;
|
|
break;
|
|
}
|
|
}
|
|
|
|
/// Returns the parametric center of the given cell shape with the given number
|
|
/// of points.
|
|
///
|
|
template<typename CellShapeTag>
|
|
VTKM_EXEC_EXPORT
|
|
vtkm::Vec<vtkm::FloatDefault,3>
|
|
ParametricCoordinatesCenter(vtkm::IdComponent numPoints,
|
|
CellShapeTag shape,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
vtkm::Vec<vtkm::FloatDefault,3> pcoords;
|
|
ParametricCoordinatesCenter(numPoints, pcoords, shape, worklet);
|
|
return pcoords;
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
template<typename ParametricCoordType>
|
|
VTKM_EXEC_EXPORT
|
|
void ParametricCoordinatesPoint(vtkm::IdComponent,
|
|
vtkm::IdComponent,
|
|
vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagEmpty,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
worklet.RaiseError("Empty cell has no points.");
|
|
pcoords[0] = pcoords[1] = pcoords[2] = 0;
|
|
}
|
|
|
|
template<typename ParametricCoordType>
|
|
VTKM_EXEC_EXPORT
|
|
void ParametricCoordinatesPoint(vtkm::IdComponent numPoints,
|
|
vtkm::IdComponent pointIndex,
|
|
vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagVertex,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
VTKM_ASSERT_EXEC(numPoints == 1, worklet);
|
|
VTKM_ASSERT_EXEC(pointIndex == 0, worklet);
|
|
pcoords[0] = 0;
|
|
pcoords[1] = 0;
|
|
pcoords[2] = 0;
|
|
}
|
|
|
|
template<typename ParametricCoordType>
|
|
VTKM_EXEC_EXPORT
|
|
void ParametricCoordinatesPoint(vtkm::IdComponent numPoints,
|
|
vtkm::IdComponent pointIndex,
|
|
vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagLine,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
VTKM_ASSERT_EXEC(numPoints == 2, worklet);
|
|
VTKM_ASSERT_EXEC((pointIndex >= 0) && (pointIndex < 2), worklet);
|
|
|
|
pcoords[0] = static_cast<ParametricCoordType>(pointIndex);
|
|
pcoords[1] = 0;
|
|
pcoords[2] = 0;
|
|
}
|
|
|
|
template<typename ParametricCoordType>
|
|
VTKM_EXEC_EXPORT
|
|
void ParametricCoordinatesPoint(vtkm::IdComponent numPoints,
|
|
vtkm::IdComponent pointIndex,
|
|
vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagTriangle,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
VTKM_ASSERT_EXEC(numPoints == 3, worklet);
|
|
|
|
switch (pointIndex)
|
|
{
|
|
case 0: pcoords[0] = 0; pcoords[1] = 0; break;
|
|
case 1: pcoords[0] = 1; pcoords[1] = 0; break;
|
|
case 2: pcoords[0] = 0; pcoords[1] = 1; break;
|
|
default: worklet.RaiseError("Bad point index.");
|
|
pcoords[0] = pcoords[1] = pcoords[2] = 0;
|
|
break;
|
|
}
|
|
pcoords[2] = 0;
|
|
}
|
|
|
|
template<typename ParametricCoordType>
|
|
VTKM_EXEC_EXPORT
|
|
void ParametricCoordinatesPoint(vtkm::IdComponent numPoints,
|
|
vtkm::IdComponent pointIndex,
|
|
vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagPolygon,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
VTKM_ASSERT_EXEC(numPoints > 0, worklet);
|
|
|
|
switch (numPoints)
|
|
{
|
|
case 1:
|
|
ParametricCoordinatesPoint(numPoints,
|
|
pointIndex,
|
|
pcoords,
|
|
vtkm::CellShapeTagVertex(),
|
|
worklet);
|
|
return;
|
|
case 2:
|
|
ParametricCoordinatesPoint(numPoints,
|
|
pointIndex,
|
|
pcoords,
|
|
vtkm::CellShapeTagLine(),
|
|
worklet);
|
|
return;
|
|
case 3:
|
|
ParametricCoordinatesPoint(numPoints,
|
|
pointIndex,
|
|
pcoords,
|
|
vtkm::CellShapeTagTriangle(),
|
|
worklet);
|
|
return;
|
|
case 4:
|
|
ParametricCoordinatesPoint(numPoints,
|
|
pointIndex,
|
|
pcoords,
|
|
vtkm::CellShapeTagQuad(),
|
|
worklet);
|
|
return;
|
|
}
|
|
|
|
// If we are here, then numPoints >= 5.
|
|
|
|
const ParametricCoordType angle =
|
|
static_cast<ParametricCoordType>(pointIndex*2*vtkm::Pi()/numPoints);
|
|
|
|
pcoords[0] = 0.5f*(vtkm::Cos(angle)+1);
|
|
pcoords[1] = 0.5f*(vtkm::Sin(angle)+1);
|
|
pcoords[2] = 0;
|
|
}
|
|
|
|
template<typename ParametricCoordType>
|
|
VTKM_EXEC_EXPORT
|
|
void ParametricCoordinatesPoint(vtkm::IdComponent numPoints,
|
|
vtkm::IdComponent pointIndex,
|
|
vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagQuad,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
VTKM_ASSERT_EXEC(numPoints == 4, worklet);
|
|
|
|
switch (pointIndex)
|
|
{
|
|
case 0: pcoords[0] = 0; pcoords[1] = 0; break;
|
|
case 1: pcoords[0] = 1; pcoords[1] = 0; break;
|
|
case 2: pcoords[0] = 1; pcoords[1] = 1; break;
|
|
case 3: pcoords[0] = 0; pcoords[1] = 1; break;
|
|
default:
|
|
worklet.RaiseError("Bad point index.");
|
|
pcoords[0] = pcoords[1] = pcoords[2] = 0;
|
|
break;
|
|
}
|
|
pcoords[2] = 0;
|
|
}
|
|
|
|
template<typename ParametricCoordType>
|
|
VTKM_EXEC_EXPORT
|
|
void ParametricCoordinatesPoint(vtkm::IdComponent numPoints,
|
|
vtkm::IdComponent pointIndex,
|
|
vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagTetra,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
VTKM_ASSERT_EXEC(numPoints == 4, worklet);
|
|
|
|
switch (pointIndex)
|
|
{
|
|
case 0: pcoords[0] = 0; pcoords[1] = 0; pcoords[2] = 0; break;
|
|
case 1: pcoords[0] = 1; pcoords[1] = 0; pcoords[2] = 0; break;
|
|
case 2: pcoords[0] = 0; pcoords[1] = 1; pcoords[2] = 0; break;
|
|
case 3: pcoords[0] = 0; pcoords[1] = 0; pcoords[2] = 1; break;
|
|
default:
|
|
worklet.RaiseError("Bad point index.");
|
|
pcoords[0] = pcoords[1] = pcoords[2] = 0;
|
|
break;
|
|
}
|
|
}
|
|
|
|
template<typename ParametricCoordType>
|
|
VTKM_EXEC_EXPORT
|
|
void ParametricCoordinatesPoint(vtkm::IdComponent numPoints,
|
|
vtkm::IdComponent pointIndex,
|
|
vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagHexahedron,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
VTKM_ASSERT_EXEC(numPoints == 8, worklet);
|
|
|
|
switch (pointIndex)
|
|
{
|
|
case 0: pcoords[0] = 0; pcoords[1] = 0; pcoords[2] = 0; break;
|
|
case 1: pcoords[0] = 1; pcoords[1] = 0; pcoords[2] = 0; break;
|
|
case 2: pcoords[0] = 1; pcoords[1] = 1; pcoords[2] = 0; break;
|
|
case 3: pcoords[0] = 0; pcoords[1] = 1; pcoords[2] = 0; break;
|
|
case 4: pcoords[0] = 0; pcoords[1] = 0; pcoords[2] = 1; break;
|
|
case 5: pcoords[0] = 1; pcoords[1] = 0; pcoords[2] = 1; break;
|
|
case 6: pcoords[0] = 1; pcoords[1] = 1; pcoords[2] = 1; break;
|
|
case 7: pcoords[0] = 0; pcoords[1] = 1; pcoords[2] = 1; break;
|
|
default:
|
|
worklet.RaiseError("Bad point index.");
|
|
pcoords[0] = pcoords[1] = pcoords[2] = 0;
|
|
break;
|
|
}
|
|
}
|
|
|
|
template<typename ParametricCoordType>
|
|
VTKM_EXEC_EXPORT
|
|
void ParametricCoordinatesPoint(vtkm::IdComponent numPoints,
|
|
vtkm::IdComponent pointIndex,
|
|
vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagWedge,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
VTKM_ASSERT_EXEC(numPoints == 6, worklet);
|
|
|
|
switch (pointIndex)
|
|
{
|
|
case 0: pcoords[0] = 0; pcoords[1] = 0; pcoords[2] = 0; break;
|
|
case 1: pcoords[0] = 0; pcoords[1] = 1; pcoords[2] = 0; break;
|
|
case 2: pcoords[0] = 1; pcoords[1] = 0; pcoords[2] = 0; break;
|
|
case 3: pcoords[0] = 0; pcoords[1] = 0; pcoords[2] = 1; break;
|
|
case 4: pcoords[0] = 0; pcoords[1] = 1; pcoords[2] = 1; break;
|
|
case 5: pcoords[0] = 1; pcoords[1] = 0; pcoords[2] = 1; break;
|
|
default:
|
|
worklet.RaiseError("Bad point index.");
|
|
pcoords[0] = pcoords[1] = pcoords[2] = 0;
|
|
break;
|
|
}
|
|
}
|
|
|
|
template<typename ParametricCoordType>
|
|
VTKM_EXEC_EXPORT
|
|
void ParametricCoordinatesPoint(vtkm::IdComponent numPoints,
|
|
vtkm::IdComponent pointIndex,
|
|
vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagPyramid,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
VTKM_ASSERT_EXEC(numPoints == 5, worklet);
|
|
|
|
switch (pointIndex)
|
|
{
|
|
case 0: pcoords[0] = 0; pcoords[1] = 0; pcoords[2] = 0; break;
|
|
case 1: pcoords[0] = 1; pcoords[1] = 0; pcoords[2] = 0; break;
|
|
case 2: pcoords[0] = 1; pcoords[1] = 1; pcoords[2] = 0; break;
|
|
case 3: pcoords[0] = 0; pcoords[1] = 1; pcoords[2] = 0; break;
|
|
case 4: pcoords[0] = 0.5; pcoords[1] = 0.5; pcoords[2] = 1; break;
|
|
default:
|
|
worklet.RaiseError("Bad point index.");
|
|
pcoords[0] = pcoords[1] = pcoords[2] = 0;
|
|
break;
|
|
}
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
/// Returns the parametric coordinate of a cell point of the given shape with
|
|
/// the given number of points.
|
|
///
|
|
template<typename ParametricCoordType>
|
|
VTKM_EXEC_EXPORT
|
|
void ParametricCoordinatesPoint(vtkm::IdComponent numPoints,
|
|
vtkm::IdComponent pointIndex,
|
|
vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagGeneric shape,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
switch (shape.Id)
|
|
{
|
|
vtkmGenericCellShapeMacro(ParametricCoordinatesPoint(numPoints,
|
|
pointIndex,
|
|
pcoords,
|
|
CellShapeTag(),
|
|
worklet));
|
|
default:
|
|
worklet.RaiseError("Bad shape given to ParametricCoordinatesPoint.");
|
|
pcoords[0] = pcoords[1] = pcoords[2] = 0;
|
|
break;
|
|
}
|
|
}
|
|
|
|
/// Returns the parametric coordinate of a cell point of the given shape with
|
|
/// the given number of points.
|
|
///
|
|
template<typename CellShapeTag>
|
|
VTKM_EXEC_EXPORT
|
|
vtkm::Vec<vtkm::FloatDefault,3>
|
|
ParametricCoordinatesPoint(vtkm::IdComponent numPoints,
|
|
vtkm::IdComponent pointIndex,
|
|
CellShapeTag shape,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
vtkm::Vec<vtkm::FloatDefault,3> pcoords;
|
|
ParametricCoordinatesPoint(numPoints, pointIndex, pcoords, shape, worklet);
|
|
return pcoords;
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
template<typename WorldCoordVector, typename PCoordType, typename CellShapeTag>
|
|
VTKM_EXEC_EXPORT
|
|
typename WorldCoordVector::ComponentType
|
|
ParametricCoordinatesToWorldCoordinates(
|
|
const WorldCoordVector &pointWCoords,
|
|
const vtkm::Vec<PCoordType,3> &pcoords,
|
|
CellShapeTag shape,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
return vtkm::exec::CellInterpolate(pointWCoords, pcoords, shape, worklet);
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
|
|
namespace detail {
|
|
|
|
template<typename WorldCoordVector, typename CellShapeTag>
|
|
class JacobianFunctorQuad
|
|
{
|
|
typedef typename WorldCoordVector::ComponentType::ComponentType T;
|
|
typedef vtkm::Vec<T,2> Vector2;
|
|
typedef vtkm::Matrix<T,2,2> Matrix2x2;
|
|
typedef vtkm::exec::internal::Space2D<T> SpaceType;
|
|
|
|
const WorldCoordVector *PointWCoords;
|
|
const SpaceType *Space;
|
|
|
|
public:
|
|
VTKM_EXEC_EXPORT
|
|
JacobianFunctorQuad(
|
|
const WorldCoordVector *pointWCoords,
|
|
const SpaceType *space)
|
|
: PointWCoords(pointWCoords), Space(space)
|
|
{ }
|
|
|
|
VTKM_EXEC_EXPORT
|
|
Matrix2x2 operator()(const Vector2 &pcoords) const
|
|
{
|
|
Matrix2x2 jacobian;
|
|
vtkm::exec::internal::JacobianFor2DCell(
|
|
*this->PointWCoords,
|
|
vtkm::Vec<T,3>(pcoords[0],pcoords[1],0),
|
|
*this->Space,
|
|
jacobian,
|
|
CellShapeTag());
|
|
return jacobian;
|
|
}
|
|
};
|
|
|
|
template<typename WorldCoordVector, typename CellShapeTag>
|
|
class CoordinatesFunctorQuad
|
|
{
|
|
typedef typename WorldCoordVector::ComponentType::ComponentType T;
|
|
typedef vtkm::Vec<T,2> Vector2;
|
|
typedef vtkm::Vec<T,3> Vector3;
|
|
typedef vtkm::exec::internal::Space2D<T> SpaceType;
|
|
|
|
const WorldCoordVector *PointWCoords;
|
|
const SpaceType *Space;
|
|
const vtkm::exec::FunctorBase *Worklet;
|
|
|
|
public:
|
|
VTKM_EXEC_EXPORT
|
|
CoordinatesFunctorQuad(
|
|
const WorldCoordVector *pointWCoords,
|
|
const SpaceType *space,
|
|
const vtkm::exec::FunctorBase *worklet)
|
|
: PointWCoords(pointWCoords), Space(space), Worklet(worklet)
|
|
{ }
|
|
|
|
VTKM_EXEC_EXPORT
|
|
Vector2 operator()(Vector2 pcoords) const {
|
|
Vector3 pcoords3D(pcoords[0], pcoords[1], 0);
|
|
Vector3 wcoords =
|
|
vtkm::exec::ParametricCoordinatesToWorldCoordinates(
|
|
*this->PointWCoords, pcoords3D, CellShapeTag(), *this->Worklet);
|
|
return this->Space->ConvertCoordToSpace(wcoords);
|
|
}
|
|
};
|
|
|
|
template<typename WorldCoordVector, typename CellShapeTag>
|
|
class JacobianFunctor3DCell
|
|
{
|
|
typedef typename WorldCoordVector::ComponentType::ComponentType T;
|
|
typedef vtkm::Vec<T,3> Vector3;
|
|
typedef vtkm::Matrix<T,3,3> Matrix3x3;
|
|
|
|
const WorldCoordVector *PointWCoords;
|
|
|
|
public:
|
|
VTKM_EXEC_EXPORT
|
|
JacobianFunctor3DCell(
|
|
const WorldCoordVector *pointWCoords)
|
|
: PointWCoords(pointWCoords)
|
|
{ }
|
|
|
|
VTKM_EXEC_EXPORT
|
|
Matrix3x3 operator()(const Vector3 &pcoords) const
|
|
{
|
|
Matrix3x3 jacobian;
|
|
vtkm::exec::internal::JacobianFor3DCell(*this->PointWCoords,
|
|
pcoords,
|
|
jacobian,
|
|
CellShapeTag());
|
|
return jacobian;
|
|
}
|
|
};
|
|
|
|
template<typename WorldCoordVector, typename CellShapeTag>
|
|
class CoordinatesFunctor3DCell
|
|
{
|
|
typedef typename WorldCoordVector::ComponentType::ComponentType T;
|
|
typedef vtkm::Vec<T,3> Vector3;
|
|
|
|
const WorldCoordVector *PointWCoords;
|
|
const vtkm::exec::FunctorBase *Worklet;
|
|
|
|
public:
|
|
VTKM_EXEC_EXPORT
|
|
CoordinatesFunctor3DCell(const WorldCoordVector *pointWCoords,
|
|
const vtkm::exec::FunctorBase *worklet)
|
|
: PointWCoords(pointWCoords), Worklet(worklet)
|
|
{ }
|
|
|
|
VTKM_EXEC_EXPORT
|
|
Vector3 operator()(Vector3 pcoords) const {
|
|
return
|
|
vtkm::exec::ParametricCoordinatesToWorldCoordinates(
|
|
*this->PointWCoords, pcoords, CellShapeTag(), *this->Worklet);
|
|
}
|
|
};
|
|
|
|
template<typename WorldCoordVector, typename CellShapeTag>
|
|
VTKM_EXEC_EXPORT
|
|
typename WorldCoordVector::ComponentType
|
|
WorldCoordinatesToParametricCoordinates3D(
|
|
const WorldCoordVector &pointWCoords,
|
|
const typename WorldCoordVector::ComponentType &wcoords,
|
|
CellShapeTag,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
return vtkm::exec::NewtonsMethod(
|
|
JacobianFunctor3DCell<WorldCoordVector,CellShapeTag>(&pointWCoords),
|
|
CoordinatesFunctor3DCell<WorldCoordVector,CellShapeTag>(&pointWCoords, &worklet),
|
|
wcoords,
|
|
typename WorldCoordVector::ComponentType(0.5f,0.5f,0.5f));
|
|
}
|
|
|
|
|
|
} // namespace detail
|
|
|
|
//-----------------------------------------------------------------------------
|
|
template<typename WorldCoordVector>
|
|
VTKM_EXEC_EXPORT
|
|
typename WorldCoordVector::ComponentType
|
|
WorldCoordinatesToParametricCoordinates(
|
|
const WorldCoordVector &pointWCoords,
|
|
const typename WorldCoordVector::ComponentType &wcoords,
|
|
vtkm::CellShapeTagGeneric shape,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
switch (shape.Id)
|
|
{
|
|
vtkmGenericCellShapeMacro(
|
|
return WorldCoordinatesToParametricCoordinates(pointWCoords,
|
|
wcoords,
|
|
CellShapeTag(),
|
|
worklet));
|
|
default:
|
|
worklet.RaiseError("Unknown cell shape sent to world 2 parametric.");
|
|
return typename WorldCoordVector::ComponentType();
|
|
}
|
|
}
|
|
|
|
template<typename WorldCoordVector>
|
|
VTKM_EXEC_EXPORT
|
|
typename WorldCoordVector::ComponentType
|
|
WorldCoordinatesToParametricCoordinates(
|
|
const WorldCoordVector &,
|
|
const typename WorldCoordVector::ComponentType &,
|
|
vtkm::CellShapeTagEmpty,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
worklet.RaiseError("Attempted to find point coordinates in empty cell.");
|
|
return typename WorldCoordVector::ComponentType();
|
|
}
|
|
|
|
template<typename WorldCoordVector>
|
|
VTKM_EXEC_EXPORT
|
|
typename WorldCoordVector::ComponentType
|
|
WorldCoordinatesToParametricCoordinates(
|
|
const WorldCoordVector &pointWCoords,
|
|
const typename WorldCoordVector::ComponentType &,
|
|
vtkm::CellShapeTagVertex,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
VTKM_ASSERT_EXEC(pointWCoords.GetNumberOfComponents() == 1, worklet);
|
|
return typename WorldCoordVector::ComponentType(0, 0, 0);
|
|
}
|
|
|
|
template<typename WorldCoordVector>
|
|
VTKM_EXEC_EXPORT
|
|
typename WorldCoordVector::ComponentType
|
|
WorldCoordinatesToParametricCoordinates(
|
|
const WorldCoordVector &pointWCoords,
|
|
const typename WorldCoordVector::ComponentType &wcoords,
|
|
vtkm::CellShapeTagLine,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
VTKM_ASSERT_EXEC(pointWCoords.GetNumberOfComponents() == 2, worklet);
|
|
|
|
// Because this is a line, there is only one vaild parametric coordinate. Let
|
|
// vec be the vector from the first point to the second point
|
|
// (pointWCoords[1] - pointWCoords[0]), which is the direction of the line.
|
|
// dot(vec,wcoords-pointWCoords[0])/mag(vec) is the orthoginal projection of
|
|
// wcoords on the line and represents the distance between the orthoginal
|
|
// projection and pointWCoords[0]. The parametric coordinate is the fraction
|
|
// of this over the length of the segment, which is mag(vec). Thus, the
|
|
// parametric coordinate is dot(vec,wcoords-pointWCoords[0])/mag(vec)^2.
|
|
|
|
typedef typename WorldCoordVector::ComponentType Vector3;
|
|
typedef typename Vector3::ComponentType T;
|
|
|
|
Vector3 vec = pointWCoords[1] - pointWCoords[0];
|
|
T numerator = vtkm::dot(vec, wcoords - pointWCoords[0]);
|
|
T denominator = vtkm::MagnitudeSquared(vec);
|
|
|
|
return Vector3(numerator/denominator, 0, 0);
|
|
}
|
|
|
|
template<typename WorldCoordVector>
|
|
VTKM_EXEC_EXPORT
|
|
typename WorldCoordVector::ComponentType
|
|
WorldCoordinatesToParametricCoordinates(
|
|
const WorldCoordVector &pointWCoords,
|
|
const typename WorldCoordVector::ComponentType &wcoords,
|
|
vtkm::CellShapeTagTriangle,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
return vtkm::exec::internal::ReverseInterpolateTriangle(
|
|
pointWCoords, wcoords, worklet);
|
|
}
|
|
|
|
template<typename WorldCoordVector>
|
|
VTKM_EXEC_EXPORT
|
|
typename WorldCoordVector::ComponentType
|
|
WorldCoordinatesToParametricCoordinates(
|
|
const WorldCoordVector &pointWCoords,
|
|
const typename WorldCoordVector::ComponentType &wcoords,
|
|
vtkm::CellShapeTagPolygon,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
const vtkm::IdComponent numPoints = pointWCoords.GetNumberOfComponents();
|
|
VTKM_ASSERT_EXEC(numPoints > 0, worklet);
|
|
switch (numPoints)
|
|
{
|
|
case 1:
|
|
return WorldCoordinatesToParametricCoordinates(
|
|
pointWCoords, wcoords, vtkm::CellShapeTagVertex(), worklet);
|
|
case 2:
|
|
return WorldCoordinatesToParametricCoordinates(
|
|
pointWCoords, wcoords, vtkm::CellShapeTagLine(), worklet);
|
|
case 3:
|
|
return WorldCoordinatesToParametricCoordinates(
|
|
pointWCoords, wcoords, vtkm::CellShapeTagTriangle(), worklet);
|
|
case 4:
|
|
return WorldCoordinatesToParametricCoordinates(
|
|
pointWCoords, wcoords, vtkm::CellShapeTagQuad(), worklet);
|
|
}
|
|
|
|
// If we are here, then there are 5 or more points on this polygon.
|
|
|
|
// Arrange the points such that they are on the circle circumscribed in the
|
|
// unit square from 0 to 1. That is, the point are on the circle centered at
|
|
// coordinate 0.5,0.5 with radius 0.5. The polygon is divided into regions
|
|
// defined by they triangle fan formed by the points around the center. This
|
|
// is C0 continuous but not necessarily C1 continuous. It is also possible to
|
|
// have a non 1 to 1 mapping between parametric coordinates world coordinates
|
|
// if the polygon is not planar or convex.
|
|
|
|
typedef typename WorldCoordVector::ComponentType WCoordType;
|
|
|
|
// Find the position of the center point.
|
|
WCoordType wcoordCenter = pointWCoords[0];
|
|
for (vtkm::IdComponent pointIndex = 1; pointIndex < numPoints; pointIndex++)
|
|
{
|
|
wcoordCenter = wcoordCenter + pointWCoords[pointIndex];
|
|
}
|
|
wcoordCenter = wcoordCenter*WCoordType(1.0f/numPoints);
|
|
|
|
// Find the normal vector to the polygon. If the polygon is planar, convex,
|
|
// and in general position, any three points will give a normal in the same
|
|
// direction. Although not perfectly robust, we can reduce the effect of
|
|
// non-planar, non-convex, or degenerate polygons by picking three points
|
|
// topologically far from each other. Note that we do not care about the
|
|
// length of the normal in this case.
|
|
WCoordType polygonNormal;
|
|
{
|
|
WCoordType vec1 = pointWCoords[numPoints/3] - pointWCoords[0];
|
|
WCoordType vec2 = pointWCoords[2*numPoints/3] - pointWCoords[1];
|
|
polygonNormal = vtkm::Cross(vec1, vec2);
|
|
}
|
|
|
|
// Find which triangle wcoords is located in. We do this by defining the
|
|
// equations for the planes through the radial edges and perpendicular to the
|
|
// polygon. The point is in the triangle if it is on the correct side of both
|
|
// planes.
|
|
vtkm::IdComponent firstPointIndex;
|
|
vtkm::IdComponent secondPointIndex;
|
|
bool foundTriangle = false;
|
|
for (firstPointIndex = 0; firstPointIndex < numPoints-1; firstPointIndex++)
|
|
{
|
|
WCoordType vecInPlane = pointWCoords[firstPointIndex] - wcoordCenter;
|
|
WCoordType planeNormal = vtkm::Cross(polygonNormal, vecInPlane);
|
|
typename WCoordType::ComponentType planeOffset =
|
|
vtkm::dot(planeNormal,wcoordCenter);
|
|
if (vtkm::dot(planeNormal,wcoords) < planeOffset)
|
|
{
|
|
// wcoords on wrong side of plane, thus outside of triangle
|
|
continue;
|
|
}
|
|
|
|
secondPointIndex = firstPointIndex+1;
|
|
vecInPlane = pointWCoords[secondPointIndex] - wcoordCenter;
|
|
planeNormal = vtkm::Cross(polygonNormal, vecInPlane);
|
|
planeOffset = vtkm::dot(planeNormal,wcoordCenter);
|
|
if (vtkm::dot(planeNormal,wcoords) > planeOffset)
|
|
{
|
|
// wcoords on wrong side of plane, thus outside of triangle
|
|
continue;
|
|
}
|
|
|
|
foundTriangle = true;
|
|
break;
|
|
}
|
|
if (!foundTriangle)
|
|
{
|
|
// wcoord was outside of all triangles we checked. It must be inside the
|
|
// one triangle we did not check (the one between the first and last
|
|
// polygon points).
|
|
firstPointIndex = numPoints-1;
|
|
secondPointIndex = 0;
|
|
}
|
|
|
|
// Build a structure containing the points of the triangle wcoords is in and
|
|
// use the triangle version of this function to find the parametric
|
|
// coordinates.
|
|
vtkm::Vec<WCoordType,3> triangleWCoords;
|
|
triangleWCoords[0] = wcoordCenter;
|
|
triangleWCoords[1] = pointWCoords[firstPointIndex];
|
|
triangleWCoords[2] = pointWCoords[secondPointIndex];
|
|
|
|
WCoordType trianglePCoords =
|
|
WorldCoordinatesToParametricCoordinates(triangleWCoords,
|
|
wcoords,
|
|
vtkm::CellShapeTagTriangle(),
|
|
worklet);
|
|
|
|
// trianglePCoords is in the triangle's parameter space rather than the
|
|
// polygon's parameter space. We can find the polygon's parameter space by
|
|
// repurposing the ParametricCoordinatesToWorldCoordinates by using the
|
|
// polygon parametric coordinates as a proxy for world coordinates.
|
|
triangleWCoords[0] = WCoordType(0.5f, 0.5f, 0);
|
|
ParametricCoordinatesPoint(numPoints,
|
|
firstPointIndex,
|
|
triangleWCoords[1],
|
|
vtkm::CellShapeTagPolygon(),
|
|
worklet);
|
|
ParametricCoordinatesPoint(numPoints,
|
|
secondPointIndex,
|
|
triangleWCoords[2],
|
|
vtkm::CellShapeTagPolygon(),
|
|
worklet);
|
|
return ParametricCoordinatesToWorldCoordinates(triangleWCoords,
|
|
trianglePCoords,
|
|
vtkm::CellShapeTagTriangle(),
|
|
worklet);
|
|
}
|
|
|
|
template<typename WorldCoordVector>
|
|
VTKM_EXEC_EXPORT
|
|
typename WorldCoordVector::ComponentType
|
|
WorldCoordinatesToParametricCoordinates(
|
|
const WorldCoordVector &pointWCoords,
|
|
const typename WorldCoordVector::ComponentType &wcoords,
|
|
vtkm::CellShapeTagQuad,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
VTKM_ASSERT_EXEC(pointWCoords.GetNumberOfComponents() == 4, worklet);
|
|
|
|
typedef typename WorldCoordVector::ComponentType::ComponentType T;
|
|
typedef vtkm::Vec<T,2> Vector2;
|
|
typedef vtkm::Vec<T,3> Vector3;
|
|
|
|
// We have an underdetermined system in 3D, so create a 2D space in the
|
|
// plane that the polygon sits.
|
|
vtkm::exec::internal::Space2D<T> space(
|
|
pointWCoords[0], pointWCoords[1], pointWCoords[3]);
|
|
|
|
Vector2 pcoords =
|
|
vtkm::exec::NewtonsMethod(
|
|
detail::JacobianFunctorQuad<WorldCoordVector,vtkm::CellShapeTagQuad>(&pointWCoords, &space),
|
|
detail::CoordinatesFunctorQuad<WorldCoordVector,vtkm::CellShapeTagQuad>(&pointWCoords, &space, &worklet),
|
|
space.ConvertCoordToSpace(wcoords),
|
|
Vector2(0.5f, 0.5f));
|
|
|
|
return Vector3(pcoords[0], pcoords[1], 0);
|
|
}
|
|
|
|
template<typename WorldCoordVector>
|
|
VTKM_EXEC_EXPORT
|
|
typename WorldCoordVector::ComponentType
|
|
WorldCoordinatesToParametricCoordinates(
|
|
const WorldCoordVector &pointWCoords,
|
|
const typename WorldCoordVector::ComponentType &wcoords,
|
|
vtkm::CellShapeTagTetra,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
VTKM_ASSERT_EXEC(pointWCoords.GetNumberOfComponents() == 4, worklet);
|
|
|
|
// We solve the world to parametric coordinates problem for tetrahedra
|
|
// similarly to that for triangles. Before understanding this code, you
|
|
// should understand the triangle code (in ReverseInterpolateTriangle in
|
|
// CellInterpolate.h). Go ahead. Read it now.
|
|
//
|
|
// The tetrahedron code is an obvious extension of the triangle code by
|
|
// considering the parallelpiped formed by wcoords and p0 of the triangle
|
|
// and the three adjacent faces. This parallelpiped is equivalent to the
|
|
// axis-aligned cuboid anchored at the origin of parametric space.
|
|
//
|
|
// Just like the triangle, we compute the parametric coordinate for each axis
|
|
// by intersecting a plane with each edge emanating from p0. The plane is
|
|
// defined by the one that goes through wcoords (duh) and is parallel to the
|
|
// plane formed by the other two edges emanating from p0 (as dictated by the
|
|
// aforementioned parallelpiped).
|
|
//
|
|
// In review, by parameterizing the line by fraction of distance the distance
|
|
// from p0 to the adjacent point (which is itself the parametric coordinate
|
|
// we are after), we get the following definition for the intersection.
|
|
//
|
|
// d = dot((wcoords - p0), planeNormal)/dot((p1-p0), planeNormal)
|
|
//
|
|
|
|
typedef typename WorldCoordVector::ComponentType Vector3;
|
|
|
|
Vector3 pcoords;
|
|
|
|
const Vector3 vec0 = pointWCoords[1] - pointWCoords[0];
|
|
const Vector3 vec1 = pointWCoords[2] - pointWCoords[0];
|
|
const Vector3 vec2 = pointWCoords[3] - pointWCoords[0];
|
|
const Vector3 coordVec = wcoords - pointWCoords[0];
|
|
|
|
Vector3 planeNormal = vtkm::Cross(vec1, vec2);
|
|
pcoords[0] = vtkm::dot(coordVec, planeNormal)/vtkm::dot(vec0, planeNormal);
|
|
|
|
planeNormal = vtkm::Cross(vec0, vec2);
|
|
pcoords[1] = vtkm::dot(coordVec, planeNormal)/vtkm::dot(vec1, planeNormal);
|
|
|
|
planeNormal = vtkm::Cross(vec0, vec1);
|
|
pcoords[2] = vtkm::dot(coordVec, planeNormal)/vtkm::dot(vec2, planeNormal);
|
|
|
|
return pcoords;
|
|
}
|
|
|
|
|
|
template<typename WorldCoordVector>
|
|
VTKM_EXEC_EXPORT
|
|
typename WorldCoordVector::ComponentType
|
|
WorldCoordinatesToParametricCoordinates(
|
|
const WorldCoordVector &pointWCoords,
|
|
const typename WorldCoordVector::ComponentType &wcoords,
|
|
vtkm::CellShapeTagHexahedron,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
VTKM_ASSERT_EXEC(pointWCoords.GetNumberOfComponents() == 8, worklet);
|
|
|
|
return detail::WorldCoordinatesToParametricCoordinates3D(
|
|
pointWCoords, wcoords, vtkm::CellShapeTagHexahedron(), worklet);
|
|
}
|
|
|
|
template<typename WorldCoordVector>
|
|
VTKM_EXEC_EXPORT
|
|
typename WorldCoordVector::ComponentType
|
|
WorldCoordinatesToParametricCoordinates(
|
|
const WorldCoordVector &pointWCoords,
|
|
const typename WorldCoordVector::ComponentType &wcoords,
|
|
vtkm::CellShapeTagWedge,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
VTKM_ASSERT_EXEC(pointWCoords.GetNumberOfComponents() == 6, worklet);
|
|
|
|
return detail::WorldCoordinatesToParametricCoordinates3D(
|
|
pointWCoords, wcoords, vtkm::CellShapeTagWedge(), worklet);
|
|
}
|
|
|
|
|
|
template<typename WorldCoordVector>
|
|
VTKM_EXEC_EXPORT
|
|
typename WorldCoordVector::ComponentType
|
|
WorldCoordinatesToParametricCoordinates(
|
|
const WorldCoordVector &pointWCoords,
|
|
const typename WorldCoordVector::ComponentType &wcoords,
|
|
vtkm::CellShapeTagPyramid,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
VTKM_ASSERT_EXEC(pointWCoords.GetNumberOfComponents() == 5, worklet);
|
|
|
|
return detail::WorldCoordinatesToParametricCoordinates3D(
|
|
pointWCoords, wcoords, vtkm::CellShapeTagPyramid(), worklet);
|
|
}
|
|
|
|
|
|
}
|
|
} // namespace vtkm::exec
|
|
|
|
#endif //vtk_m_exec_ParametricCoordinates_h
|