Commit e7a0931f authored by Steven Walton's avatar Steven Walton

Stretch Done

Used autos, should probably use something different but it should work
since we are using 11
parent 0261c8bc
......@@ -39,6 +39,9 @@
#include "vtkm/VectorAnalysis.h"
#include "vtkm/exec/FunctorBase.h"
static constexpr FloatType FLOAT_MAX = vtkm::Infinity<FloatType>();
static constexpr FloatType FLOAT_MIN = vtkm::NegativeInfinity<FloatType>();
template <typename OutType, typename PointCoordVecType, typename CellShapeType>
VTKM_EXEC OutType CellStretchMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
......@@ -46,7 +49,7 @@ VTKM_EXEC OutType CellStretchMetric(const vtkm::IdComponent& numPts,
const vtkm::exec::FunctorBase& worklet)
{
worklet.RaiseError("Shape type template must be specified to compute stretch")
return OutType(0.0);
return OutType(-1.0);
}
template<typename OutType, typename PointCoordVecType>
......@@ -55,21 +58,17 @@ VTKM_EXEC OutType CellStretchMetric(const vtkm::IdComponent& numPts,
vtlm::CellShapeTagQuad,
const vtkm::exec::FunctorBase& worklet)
{
double min_L = vtkm::Pow(pts[1][0] - pts[0][0],2) + vtkm::Pow(pts[1][1] - pts[0][1],2);
auto L0 = vtkm::MagnitudeSquared(pts[1] - pts[0]);
auto L1 = vtkm::MagnitudeSquared(pts[2] - pts[1]);
auto L2 = vtkm::MagnitudeSquared(pts[3] - pts[2]);
auto L3 = vtkm::MagnitudeSquared(pts[0] - pts[3]);
// Find the minimum length (use square of values to speed up)
for(size_t i = 1; i < numPts; ++i)
{
double L = vtkm::Pow(pts[(i+1)%numPts][0] - pts[i][0],2) + vtkm::Pow(pts[(i+1)%numPts][1] - pts[i][1],2);
if( L < min_L)
min_L = L;
}
double D_max = vtkm::Pow(pts[2][0] - pts[0][0],2) + vtkm::Pow(pts[2][1] - pts[0][1],2);
double D2 = vtkm::Pow(pts[3][0] - pts[1][0],2) + vtkm::Pow(pts[3][1] - pts[1][1],2);
if( D_max < D2 )
D_max = D2;
auto D0 = pts[2] - pts[0];
auto D1 = pts[3] - pts[1];
auto D_max = vtkm::Max(D0,D1);
if(D_max < FLOAT_MIN) return FLOAT_MAX;
q = vtkm::Sqrt(2)/D_max * vtkm::Sqrt(min_L);
return q;
return vtkm::Sqrt(2)/D_max * vtkm::Sqrt(L_min);
}
template<typename OutType, typename PointCoordVecType>
......@@ -78,34 +77,25 @@ VTKM_EXEC OutType CellStretchMetric(const vtkm::IdComponent& numPts,
vtlm::CellShapeTagHex,
const vtkm::exec::FunctorBase& worklet)
{
double min_L = vtkm::Pow(pts[1][0] - pts[0][0],2) + vtkm::Pow(pts[1][1] - pts[0][1],2);
// Find the minimum length (use square of values to speed up)
double L = 0;
for(size_t i = 1; i < numPts; ++i)
{
L = vtkm::Pow(pts[(i+1)%numPts/2][0] - pts[i][0],2) + vtkm::Pow(pts[(i+1)%numPts/2][1] - pts[i][1],2);
if( L < min_L)
min_L = L;
}
L = vtkm::Pow(pts[4][0] - pts[0][0],2) + vtkm::Pow(pts[4][1] - pts[0][1],2);
if(L < min_L)
min_L = L;
for(size_t i = numPts; i < numPts*2; ++i)
{
double L = vtkm::Pow(pts[(i+1)%numPts/2][0] - pts[i][0],2) + vtkm::Pow(pts[(i+1)%numPts/2][1] - pts[i][1],2);
if(L < min_L)
min_L = L;
}
double D0 = vtkm::Pow(pts[6][0] - pts[0][0],2) + vtkm::Pow(pts[6][1] - pts[0][1]);
double D1 = vtkm::Pow(pts[7][0] - pts[1][0],2) + vtkm::Pow(pts[7][1] - pts[1][1]);
double D2 = vtkm::Pow(pts[4][0] - pts[2][0],2) + vtkm::Pow(pts[4][1] - pts[2][1]);
double D3 = vtkm::Pow(pts[5][0] - pts[3][0],3) + vtkm::Pow(pts[5][1] - pts[3][1]);
double D_max = vtkm::Max(D0,D1,D2,D3);
q = vtkm::Sqrt(3) / D_max * vtkm::Sqrt(min_L);
return q;
auto L0 = vtkm::MagnitudeSquared(pts[1] - pts[0]);
auto L1 = vtkm::MagnitudeSquared(pts[2] - pts[1]);
auto L2 = vtkm::MagnitudeSquared(pts[3] - pts[2]);
auto L3 = vtkm::MagnitudeSquared(pts[3] - pts[0]);
auto L4 = vtkm::MagnitudeSquared(pts[4] - pts[0]);
auto L5 = vtkm::MagnitudeSquared(pts[5] - pts[1]);
auto L6 = vtkm::MagnitudeSquared(pts[6] - pts[2]);
auto L7 = vtkm::MagnitudeSquared(pts[3] - pts[3]);
auto L8 = vtkm::MagnitudeSquared(pts[5] - pts[4]);
auto L9 = vtkm::MagnitudeSquared(pts[6] - pts[5]);
auto L10 = vtkm::MagnitudeSquared(pts[7] - pts[6]);
auto L11 = vtkm::MagnitudeSquared(pts[7] - pts[4]);
auto D0 = pts[6] - pts[0];
auto D1 = pts[7] - pts[1];
auto D2 = pts[4] - pts[2];
auto D3 = pts[5] - pts[3];
auto D_max = vtkm::Max(D0,D1,D2,D3);
if(D_max < FLOAT_MIN) return FLOAT_MAX;
return vtkm::Sqrt(3) * vtkm::Sqrt(vtkm::Min(L0,L1,L2,L3,L4,L5,L6,L7,L8,L9,L10,L11)) / D_max;
}
......
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