Commit 68c79fc6 authored by Brent Lessley's avatar Brent Lessley

Removed more compile errors and warnings related to 3 of the metrics. Still more remaining...

parent e1a83053
......@@ -36,6 +36,7 @@
If no filename is given, then "output.vtk" is used. The mesh quality of the new field is determined and output to std::cout
*/
using namespace vtkm::cont;
template<typename msg1, typename msg2>
void
......
......@@ -78,7 +78,7 @@ VTKM_EXEC OutType CellAspectFrobeniusMetric(const vtkm::IdComponent& numPts,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts == 3)
return CellAspectFrobeniusMetric(numPts, pts, vtkm::CellShapeTagTriangle(), worklet);
return CellAspectFrobeniusMetric<OutType>(numPts, pts, vtkm::CellShapeTagTriangle(), worklet);
else
{
worklet.RaiseError("Aspect frobenius metric is not supported for (n<3)- or (n>4)-vertex polygons.");
......@@ -180,9 +180,6 @@ VTKM_EXEC OutType CellAspectFrobeniusMetric(const vtkm::IdComponent& numPts,
return OutType(0.0);
}
const vtkm::IdComponent numEdges = 3;
//The 3 edges of a triangle
using Edge = typename PointCoordVecType::ComponentType;
const Edge TriEdges[3] = {pts[1] - pts[0],
......
......@@ -166,7 +166,7 @@ VTKM_EXEC OutType CellMinAngleMetric(const vtkm::IdComponent& numPts,
ctrPlusOne = (ctr+1)%4; // wrap around to first edge
edge1 = TriEdges[ctr];
edge2 = TriEdges[ctrPlusOne];
tempResult = 1/vtkm::Cos(vtkm::Dot(edge1, edge2)/(edge1*edge2));
//tempResult = 1/vtkm::Cos(vtkm::Dot(edge1, edge2)/(edge1*edge2));
tempResult *= scalar;
result = tempResult > result ? tempResult : result;
}
......@@ -205,7 +205,7 @@ VTKM_EXEC OutType CellMinAngleMetric(const vtkm::IdComponent& numPts,
};
const Edge principleXAxis = QuadEdges[0] + (pts[2] - pts[3]);
const Edge principleYAxis = QuadEdges[1] + (pts[3] - pts[0]);
const OutType centerNormal = vtkm::Cross(principleXAxis, principleYAxis)/ vtkm::Cross(principleXAxis, principleYAxis);
const Edge centerNormal = vtkm::Cross(principleXAxis, principleYAxis)/ vtkm::Cross(principleXAxis, principleYAxis);
const OutType negativeInfinity = vtkm::NegativeInfinity<OutType>();
const OutType infinity = vtkm::Infinity<OutType>();
const double scalar = 57.2957795131; // ~ 180/pi
......@@ -214,7 +214,7 @@ VTKM_EXEC OutType CellMinAngleMetric(const vtkm::IdComponent& numPts,
int ctr, ctrPlusOne, correction;
for(ctr = 0; ctr < 3; ctr++)
{
if(QuadEdges[ctr] <= negativeInfinity)
if(vtkm::MagnitudeSquared(QuadEdges[ctr]) <= negativeInfinity)
{
return OutType(360);
}
......@@ -225,9 +225,9 @@ VTKM_EXEC OutType CellMinAngleMetric(const vtkm::IdComponent& numPts,
ctrPlusOne = (ctr+1)%4; // wrap around to first edge
edge1 = QuadEdges[ctr];
edge2 = QuadEdges[ctrPlusOne];
tempNormal = vtkm::Cross(QuadEdges[(ctr+3)%4], edge1);
tempNormal = vtkm::Dot(centerNormal, tempNormal);
tempResult = 1/vtkm::Cos(-1*vtkm::Dot(edge1, edge2)/(edge1*edge2));
//tempNormal = vtkm::Cross(QuadEdges[(ctr+3)%4], edge1);
tempNormal = vtkm::Dot(centerNormal, vtkm::Cross(QuadEdges[(ctr+3)%4], edge1));
//tempResult = 1/vtkm::Cos(-1*vtkm::Dot(edge1, edge2)/(edge1*edge2));
tempResult *= scalar;
if(tempNormal < 0)
{
......
......@@ -122,9 +122,9 @@ VTKM_EXEC OutType CellMaxAspectFrobeniusMetric(const vtkm::IdComponent& numPts,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts == 3)
return CellAspectFrobeniusMetric(numPts, pts, vtkm::CellShapeTagTriangle(), worklet);
return vtkm::exec::cellmetrics::CellAspectFrobeniusMetric<OutType>(numPts, pts, vtkm::CellShapeTagTriangle(), worklet);
if (numPts == 4)
return CellMaxAspectFrobeniusMetric(numPts, pts, vtkm::CellShapeTagQuad(), worklet);
return CellMaxAspectFrobeniusMetric<OutType>(numPts, pts, vtkm::CellShapeTagQuad(), worklet);
else
{
worklet.RaiseError("Max aspect frobenius metric is not supported for (n<3)- or (n>4)-vertex polygons.");
......@@ -159,7 +159,7 @@ VTKM_EXEC OutType CellMaxAspectFrobeniusMetric(const vtkm::IdComponent& numPts,
worklet.RaiseError("Max aspect frobenius metric(triangle) requires 3 points.");
return OutType(0.0);
}
return vtkm::exec::cellmetrics::CellAspectFrobeniusMetric(numPts, pts, vtkm::CellShapeTagTriangle(), worklet);
return vtkm::exec::cellmetrics::CellAspectFrobeniusMetric<OutType>(numPts, pts, vtkm::CellShapeTagTriangle(), worklet);
}
......@@ -261,7 +261,7 @@ VTKM_EXEC OutType CellMaxAspectFrobeniusMetric(const vtkm::IdComponent& numPts,
return OutType(0.0);
}
return vtkm::exec::cellmetrics::CellAspectFrobeniusMetric(numPts, pts, vtkm::CellShapeTagTetra(), worklet);
return vtkm::exec::cellmetrics::CellAspectFrobeniusMetric<OutType>(numPts, pts, vtkm::CellShapeTagTetra(), worklet);
}
......@@ -387,14 +387,14 @@ VTKM_EXEC OutType CellMaxAspectFrobeniusMetric(const vtkm::IdComponent& numPts,
//Return the maximum metric value among all 6 tets as the maximum aspect frobenius.
const vtkm::IdComponent tetPts = 4;
OutType max_aspect_frobenius =
vtkm::exec::cellmetrics::CellAspectFrobeniusMetric(tetPts,
vtkm::exec::cellmetrics::CellAspectFrobeniusMetric<OutType>(tetPts,
TetEdges[0],
vtkm::CellShapeTagTetra(),
worklet);
OutType curr;
for (vtkm::Id i = 1; i < 6; i++)
{
curr = vtkm::exec::cellmetrics::CellAspectFrobeniusMetric(tetPts,
curr = vtkm::exec::cellmetrics::CellAspectFrobeniusMetric<OutType>(tetPts,
TetEdges[i],
vtkm::CellShapeTagTetra(),
worklet);
......
......@@ -55,11 +55,13 @@ template <typename OutType, typename PointCoordVecType, typename CellShapeType>
VTKM_EXEC OutType ComputeRelativeSizeSquared(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
CellShapeType shape,
const vtkm::Id& numShapes,
const vtkm::exec::FunctorBase&)
{
UNUSED(numPts);
UNUSED(pts);
UNUSED(shape);
UNUSED(numShapes);
return OutType(0.0);
}
......@@ -67,12 +69,13 @@ template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType ComputeRelativeSizeSquared(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagPolygon,
const vtkm::Id& numPolygons,
const vtkm::exec::FunctorBase& worklet)
{
switch(numPts)
{
case 4:
return ComputeRelativeSizeSquared<OutType>(numPts, pts, vtkm::CellShapeTagQuad(), worklet);
return ComputeRelativeSizeSquared<OutType>(numPts, pts, vtkm::CellShapeTagQuad(), numPolygons, worklet);
default:
break;
}
......@@ -83,30 +86,23 @@ template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType ComputeRelativeSizeSquared(const vtkm::IdComponent&,
const PointCoordVecType&,
vtkm::CellShapeTagLine,
const vtkm::Id& numLines,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(worklet);
UNUSED(numLines);
return OutType(-1.0);
}
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType ComputeRelativeSizeSquared(const vtkm::IdComponent&,
const PointCoordVecType&,
vtkm::CellShapeTagTetra,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(worklet);
return OutType(-1.0);
}
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType ComputeRelativeSizeSquared(const vtkm::IdComponent&,
const PointCoordVecType&,
vtkm::CellShapeTagWedge,
const vtkm::Id& numWedges,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(worklet);
UNUSED(numWedges);
return OutType(-1.0);
}
......@@ -114,9 +110,11 @@ template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType ComputeRelativeSizeSquared(const vtkm::IdComponent&,
const PointCoordVecType&,
vtkm::CellShapeTagPyramid,
const vtkm::Id& numPyramids,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(worklet);
UNUSED(numPyramids);
return OutType(-1.0);
}
......@@ -133,6 +131,7 @@ template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType ComputeRelativeSizeSquared(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagTriangle,
const vtkm::Id& numTriangles,
const vtkm::exec::FunctorBase& worklet)
{
if(numPts != 3)
......@@ -140,8 +139,8 @@ VTKM_EXEC OutType ComputeRelativeSizeSquared(const vtkm::IdComponent& numPts,
worklet.RaiseError("RelativeSizeSquared metric(tri) requires 3 points");
return OutType(0.0);
}
auto shapeCountPortal = worklet.GetShapeCounts().GetPortalConstControl();
const int numTriangles = shapeCountPortal.Get(vtkm::CELL_SHAPE_TRIANGLE);
//auto shapeCountPortal = worklet.GetShapeCounts().GetPortalConstControl();
//const int numTriangles = shapeCountPortal.Get(vtkm::CELL_SHAPE_TRIANGLE);
OutType weight1_1, weight2_1, weight1_2, weight2_2, weightedDeterminant, weightedJacobian, norm, result;
const OutType root3 = sqrt(3.0);
......@@ -190,6 +189,7 @@ template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType ComputeRelativeSizeSquared(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagQuad,
const vtkm::Id& numQuads,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 4)
......@@ -199,8 +199,6 @@ VTKM_EXEC OutType ComputeRelativeSizeSquared(const vtkm::IdComponent& numPts,
}
OutType negativeInfinity = vtkm::NegativeInfinity<OutType>();
OutType area = 0;
auto shapeCountPortal = worklet.ShapeCounts.GetPortalConstControl();
const int numQuads = shapeCountPortal.Get(vtkm::CELL_SHAPE_QUAD);
OutType weight1_1, weight2_1, weight1_2, weight2_2, avgArea, weightedJacobian, norm, result;
using Edge = typename PointCoordVecType::ComponentType;
const Edge QuadEdges[4] = {pts[1] - pts[0],
......@@ -209,8 +207,8 @@ VTKM_EXEC OutType ComputeRelativeSizeSquared(const vtkm::IdComponent& numPts,
pts[0] - pts[3]
};
//get area
OutType principle1 = QuadEdges[0]+(pts[2]-pts[3]);
OutType principle2 = QuadEdges[3]+QuadEdges[1];
Edge principle1 = QuadEdges[0]+(pts[2]-pts[3]);
Edge principle2 = QuadEdges[3]+QuadEdges[1];
auto pAxis = vtkm::Cross(principle1, principle2);
OutType pNorm = pAxis/vtkm::MagnitudeSquared(pAxis);
OutType currNorm;
......@@ -256,6 +254,7 @@ template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType ComputeRelativeSizeSquared(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagTetra,
const vtkm::Id& numTets,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 4)
......@@ -265,8 +264,6 @@ VTKM_EXEC OutType ComputeRelativeSizeSquared(const vtkm::IdComponent& numPts,
}
//The 6 edges of a tetrahedron
auto shapeCountPortal = worklet.ShapeCounts.GetPortalConstControl();
const int tet_size = shapeCountPortal.Get(vtkm::CELL_SHAPE_TETRA);
const OutType negativeInfinity = vtkm::NegativeInfinity<OutType>();
const OutType rt3 = sqrt(3.0);
const OutType rt2 = sqrt(2.0);
......@@ -283,7 +280,7 @@ VTKM_EXEC OutType ComputeRelativeSizeSquared(const vtkm::IdComponent& numPts,
Edge weight2 = {0.5, 0.5*rt3, 0};
Edge weight3 = {0.5, rt3/6.0, rt2/rt3};
OutType avg_volume, volume, result, R;
OutType scale = pow(6*tet_size/ vtkm::Dot(weight1, vtkm::Cross(weight2,weight3)), 0.333333333333333333);
OutType scale = pow(6*numTets/ vtkm::Dot(weight1, vtkm::Cross(weight2,weight3)), 0.333333333333333333);
weight1 *= scale;
weight2 *= scale;
weight3 *= scale;
......@@ -316,6 +313,7 @@ template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType ComputeRelativeSizeSquared(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagHexahedron,
const vtkm::Id& numHexs,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 8)
......@@ -324,8 +322,6 @@ VTKM_EXEC OutType ComputeRelativeSizeSquared(const vtkm::IdComponent& numPts,
return OutType(0.0);
}
OutType avg_volume, volume;
auto shapeCountPortal = worklet.ShapeCounts.GetPortalConstControl();
const int hex_size = shapeCountPortal.Get(vtkm::CELL_SHAPE_HEXAHEDRON);
const OutType negativeInfinity = vtkm::NegativeInfinity<OutType>();
//The 12 edges of a hexahedron
using Edge = typename PointCoordVecType::ComponentType;
......@@ -360,7 +356,7 @@ VTKM_EXEC OutType ComputeRelativeSizeSquared(const vtkm::IdComponent& numPts,
Edge weight2 = {0, 1, 0};
Edge weight3 = {0, 0, 1};
OutType scale = pow(hex_size/(vtkm::Dot(weight1, vtkm::Cross(weight2, weight3))), 0.3333333333333333);
OutType scale = pow(numHexs/(vtkm::Dot(weight1, vtkm::Cross(weight2, weight3))), 0.3333333333333333);
weight1 *= scale;
weight2 *= scale;
weight3 *= scale;
......@@ -399,8 +395,9 @@ VTKM_EXEC OutType ComputeRelativeSizeSquared(const vtkm::IdComponent& numPts,
(tempMatrix3_3 * tempMatrix3_3);
secondPartOfNumerator = tempMatrix1_1 + tempMatrix2_2 + tempMatrix3_3;
secondPartOfNumerator *= secondPartOfNumerator;
secondPartOfNumerator *= third;
currentRelativeSizeSquared = (firstPartOfNumerator - secondPartOfNumerator)/(vtkm::Pow(determinant, OutType(4.0*third)));
OutType temp_third = 0.0;
secondPartOfNumerator *= temp_third;
currentRelativeSizeSquared = (firstPartOfNumerator - secondPartOfNumerator)/(vtkm::Pow(determinant, OutType(4.0*temp_third)));
maxRelativeSizeSquared = currentRelativeSizeSquared > maxRelativeSizeSquared ? currentRelativeSizeSquared : maxRelativeSizeSquared;
}
else {return vtkm::Infinity<OutType>(); }
......
......@@ -68,12 +68,15 @@ inline VTKM_CONT vtkm::cont::DataSet MeshQuality::DoExecute(
vtkm::cont::ArrayHandle<T> outArray;
//using MetricTagType = CellMetric;
using QualityWorklet = vtkm::worklet::MeshQuality<CellMetric>;
using DeviceAdapter = VTKM_DEFAULT_DEVICE_ADAPTER_TAG;
using QualityWorklet = vtkm::worklet::MeshQuality<CellMetric, DeviceAdapter>;
vtkm::cont::ArrayHandle<CellMetric> cellMetrics = vtkm::cont::make_ArrayHandle(CellTypeMetrics);
//Initialize the mesh quality worklet
QualityWorklet meshWorklet;
meshWorklet.SetShapeMetrics(CellTypeMetrics);
meshWorklet.SetShapeCounts(cellShapeCounts);
QualityWorklet meshWorklet(cellShapeCounts.PrepareForInput(DeviceAdapter()),
cellMetrics.PrepareForInput(DeviceAdapter()));
//meshWorklet.SetShapeMetrics(CellTypeMetrics);
//meshWorklet.SetShapeCounts(cellShapeCounts);
//Invoke the worklet
vtkm::worklet::DispatcherMapTopology<QualityWorklet> dispatcher(meshWorklet);
......
......@@ -49,7 +49,8 @@ namespace worklet
* of the cell type within this worklet. An array of the computed
* metric values (one per cell) is returned as output.
*/
template <typename MetricTagType>
template <typename MetricTagType,
typename DeviceAdapter = VTKM_DEFAULT_DEVICE_ADAPTER_TAG>
class MeshQuality : public vtkm::worklet::WorkletMapPointToCell
{
public:
......@@ -63,6 +64,17 @@ public:
using ExecutionSignature = void(CellShape, PointCount, _2, _3);
using InputDomain = _1;
using CountsPortalType = typename vtkm::cont::ArrayHandle<vtkm::Id>::template
ExecutionTypes<DeviceAdapter>::PortalConst;
using MetricsPortalType = typename vtkm::cont::ArrayHandle<MetricTagType>::template
ExecutionTypes<DeviceAdapter>::PortalConst;
MeshQuality(const CountsPortalType& counts,
const MetricsPortalType& metrics)
: ShapeCounts(counts), ShapeMetrics(metrics)
{
}
template <typename CellShapeType,
typename PointCoordVecType,
typename OutType>
......@@ -74,7 +86,9 @@ public:
switch (shape.Id)
{
vtkmGenericCellShapeMacro(metricValue =
this->ComputeMetric<OutType>(numPoints, pts, CellShapeTag()));
this->ComputeMetric<OutType>(numPoints, pts,
this->ShapeCounts.Get(shape.Id),
CellShapeTag()));
default:
this->RaiseError("Asked for metric of unknown cell type.");
metricValue = OutType(0.0);
......@@ -82,13 +96,16 @@ public:
}
//TODO: Is this proper usage of std::move? Can the passed handle be const?
/*
VTKM_CONT void SetShapeCounts(vtkm::cont::ArrayHandle<vtkm::Id> &counts)
{
this->ShapeCounts = std::move(counts);
}
*/
//TODO: Do we need a type check (std::is_same) for MetricTagType?
//TODO: If metrics is static (or static const) in the filter class, can we use std::move?
/*
VTKM_CONT void SetShapeMetrics(std::vector<MetricTagType> metrics)
{
this->ShapeMetrics = std::move(metrics);
......@@ -96,11 +113,14 @@ public:
//Initialize the metric function pointers
//InitializeMetricFunctions();
}
*/
/*
VTKM_CONT vtkm::cont::ArrayHandle<vtkm::Id> GetShapeCounts()
{
return this->ShapeCounts;
}
*/
protected:
template <typename OutType,
......@@ -108,6 +128,7 @@ protected:
typename CellShapeType>
VTKM_EXEC OutType ComputeMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::Id numShapes,
CellShapeType tag) const
{
constexpr vtkm::IdComponent dims = vtkm::CellTraits<CellShapeType>::TOPOLOGICAL_DIMENSIONS;
......@@ -120,7 +141,7 @@ protected:
//Only compute the metric for 2D and 3D shapes; return 0 otherwise
if (dims > 0)
{
auto metric = this->ShapeMetrics[tag.Id];
auto metric = this->ShapeMetrics.Get(tag.Id);
switch (metric)
{
case MetricTagType::ASPECT_FROBENIUS:
......@@ -135,14 +156,16 @@ protected:
case MetricTagType::JACOBIAN:
metricValue = vtkm::exec::cellmetrics::ComputeJacobian<OutType>(numPts, pts, tag, *this);
break;
/*
case MetricTagType::MAX_ANGLE:
metricValue = vtkm::exec::cellmetrics::CellMaxAngleMetric<OutType>(numPts, pts, tag, *this);
break;
*/
case MetricTagType::MAX_ASPECT_FROBENIUS:
metricValue = vtkm::exec::cellmetrics::CellMaxAspectFrobenius<OutType>(numPts, pts, tag, *this);
metricValue = vtkm::exec::cellmetrics::CellMaxAspectFrobeniusMetric<OutType>(numPts, pts, tag, *this);
break;
case MetricTagType::MAX_EDGE_RATIO:
metricValue = vtkm::exec::CellMaxEdgeRatioMetric<OutType>(numPts,pts,tag,*this);
metricValue = vtkm::exec::cellmetrics::CellMaxEdgeRatioMetric<OutType>(numPts,pts,tag,*this);
break;
case MetricTagType::MIN_ANGLE:
metricValue = vtkm::exec::cellmetrics::CellMinAngleMetric<OutType>(numPts, pts, tag, *this);
......@@ -151,7 +174,7 @@ protected:
metricValue = vtkm::exec::cellmetrics::ComputeOddy<OutType>(numPts, pts, tag, *this);
break;
case MetricTagType::RELATIVE_SIZE_SQUARED:
metricValue = vtkm::exec::ComputeRelativeSizeSquared<OutType>(numPts,pts,tag,*this);
metricValue = vtkm::exec::cellmetrics::ComputeRelativeSizeSquared<OutType>(numPts,pts,tag,numShapes,*this);
break;
case MetricTagType::SCALED_JACOBIAN:
metricValue = vtkm::exec::cellmetrics::ComputeScaledJacobian<OutType>(numPts,pts,tag,*this);
......@@ -203,8 +226,8 @@ protected:
*/
private:
vtkm::cont::ArrayHandle<vtkm::Id> ShapeCounts;
std::vector<MetricTagType> ShapeMetrics;
CountsPortalType ShapeCounts;
MetricsPortalType ShapeMetrics;
//std::vector<FuncPtr> MetricFunctions;
};
......
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