From 880d8a989e50441e451f0cc354f0b49bb2c5efc5 Mon Sep 17 00:00:00 2001 From: David Thompson Date: Fri, 8 Jun 2018 00:27:30 -0400 Subject: [PATCH] Add `vtkm/Geometry.h` and test it. This commit adds several geometric constructs to vtk-m in the `vtkm/Geometry.h` header. They may be used from both the execution and control environments. We also add methods to perform projection and Gram-Schmidt orthonormalization to `vtkm/VectorAnalysis.h`. See `docs/changelog/geometry.md` included in this commit for more information. --- docs/changelog/geometry.md | 75 ++ vtkm/CMakeLists.txt | 14 +- vtkm/Geometry.h | 425 ++++++++++++ vtkm/Geometry.hxx | 647 ++++++++++++++++++ vtkm/Math.h | 12 + vtkm/Math.h.in | 12 + vtkm/TypeListTag.h | 13 + vtkm/VectorAnalysis.h | 78 +++ vtkm/cont/cuda/testing/CMakeLists.txt | 1 + .../cont/cuda/testing/UnitTestCudaGeometry.cu | 38 + vtkm/cont/serial/testing/CMakeLists.txt | 1 + .../serial/testing/UnitTestSerialGeometry.cxx | 38 + vtkm/testing/CMakeLists.txt | 1 + vtkm/testing/TestingGeometry.h | 498 ++++++++++++++ vtkm/testing/UnitTestVectorAnalysis.cxx | 48 ++ 15 files changed, 1899 insertions(+), 2 deletions(-) create mode 100644 docs/changelog/geometry.md create mode 100644 vtkm/Geometry.h create mode 100644 vtkm/Geometry.hxx create mode 100644 vtkm/cont/cuda/testing/UnitTestCudaGeometry.cu create mode 100644 vtkm/cont/serial/testing/UnitTestSerialGeometry.cxx create mode 100644 vtkm/testing/TestingGeometry.h diff --git a/docs/changelog/geometry.md b/docs/changelog/geometry.md new file mode 100644 index 000000000..075328688 --- /dev/null +++ b/docs/changelog/geometry.md @@ -0,0 +1,75 @@ +# New geometry classes and header. + +There are now some additional structures available both +the control and execution environments for representing +geometric entities (mostly of dimensions 2 and 3). +These new structures are now in `vtkm/Geometry.h` and +demonstrated/tested in `vtkm/testing/TestingGeometry.h`: + ++ `Ray`. + Instances of this struct represent a semi-infinite line + segment in a 2-D plane or in a 3-D space, depending on the + integer dimension specified as a template parameter. + Its state is the point at the start of the ray (`Origin`) + plus the ray's `Direction`, a unit-length vector. + If the third template parameter (IsTwoSided) is true, then + the ray serves as an infinite line. Otherwise, the ray will + only report intersections in its positive halfspace. ++ `LineSegment`. + Instances of this struct represent a finite line segment + in a 2-D plane or in a 3-D space, depending on the integer + dimension specified as a template parameter. + Its state is the coordinates of its `Endpoints`. ++ `Plane`. + Instances of this struct represent a plane in 3-D. + Its state is the coordinates of a base point (`Origin`) and + a unit-length normal vector (`Normal`). ++ `Sphere`. + Instances of this struct represent a *d*-dimensional sphere. + Its state is the coordinates of its center plus a radius. + It is also aliased with a `using` statment to `Circle` + for the specific case of 2-D. + +These structures provide useful queries and generally +interact with one another. +For instance, it is possible to intersect lines and planes +and compute distances. + +For ease of use, there are also several `using` statements +that alias these geometric structures to names that specialize +them for a particular dimension or other template parameter. +As an example, `Ray` is aliased +to `Line` and `Ray` +is aliased to `Line3` and `Ray` +is aliased to `Line3d`. + +## Design patterns + +If you plan to add a new geometric entity type, +please adopt these conventions: + ++ Each geometric entity may be default-constructed. + The default constructor will initialize the state to some + valid unit-length entity, usually with some part of + its state at the origin of the coordinate system. ++ Entities may always be constructed by passing in values + for their internal state. + Alternate construction methods are declared as free functions + such as `make_CircleFrom3Points()` ++ Use template metaprogramming to make methods available + only when the template dimension gives them semantic meaning. + For example, a 2-D line segment's perpendicular bisector + is another line segment, but a 3-D line segment's perpendicular + line segment is a plane. + Note how this is accomplished and apply this pattern to + new geometric entities or new methods on existing entities. ++ Some entities may have invalid state. + If this is possible, the entity will have an `IsValid()` method. + For example, a sphere may be invalid because the user or some + construction technique specified a zero or negative radius. ++ When signed distance is semantically meaningful, provide it + in favor of or in addition to unsigned distance. ++ Accept a tolerance parameter when appropriate, + but provide a sensible default value. + You may want to perform exact arithmetic versions of tests, + but please provide fast, tolerance-based versions as well. diff --git a/vtkm/CMakeLists.txt b/vtkm/CMakeLists.txt index be9ff8b80..da87c186d 100644 --- a/vtkm/CMakeLists.txt +++ b/vtkm/CMakeLists.txt @@ -35,6 +35,7 @@ set(headers CellShape.h CellTraits.h Flags.h + Geometry.h Hash.h ImplicitFunction.h ListTag.h @@ -60,11 +61,20 @@ set(headers VecVariable.h VirtualObjectBase.h UnaryPredicates.h - ) +) + +set(template_sources + Geometry.hxx +) vtkm_pyexpander_generated_file(Math.h) -vtkm_declare_headers(${headers}) +vtkm_declare_headers( + ${headers} + ${template_sources} + EXCLUDE_FROM_TESTING + ${template_sources} +) #----------------------------------------------------------------------------- #first add all the components vtkm that are shared between control and exec diff --git a/vtkm/Geometry.h b/vtkm/Geometry.h new file mode 100644 index 000000000..df4286c38 --- /dev/null +++ b/vtkm/Geometry.h @@ -0,0 +1,425 @@ +//============================================================================ +// 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 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +// Copyright 2014 UT-Battelle, LLC. +// Copyright 2014 Los Alamos National Security. +// +// Under the terms of Contract DE-NA0003525 with NTESS, +// 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_Geometry_h +#define vtk_m_Geometry_h + +#include + +namespace vtkm +{ + +// Forward declarations of geometric types: +template +struct Ray; +template +struct LineSegment; +template +struct Plane; +template +struct Sphere; + +/// Represent an infinite or semi-infinite line segment with a point and a direction +/// +/// The \a IsTwoSided template parameter indicates whether the class represents +/// an infinite line extending in both directions from the base point or a +/// semi-infinite ray extending only in the positive direction from the base points. +template +struct Ray +{ + static constexpr int Dimension = Dim; + using Vector = vtkm::Vec; + static constexpr bool TwoSided = IsTwoSided; + + Vector Origin; + Vector Direction; // always stored as a unit length. + + /// Construct a default 2-D ray, from (0,0) pointing along the +x axis. + template ::type = 0> + VTKM_EXEC_CONT Ray(); + + /// Construct a default 3-D ray from (0,0,0) pointing along the +x axis. + template ::type = 0> + VTKM_EXEC_CONT Ray(); + + /// Construct a ray from a point and direction. + VTKM_EXEC_CONT + Ray(const Vector& point, const Vector& direction); + + /// Construct a ray from a line segment. + VTKM_EXEC_CONT + Ray(const LineSegment& segment); + + /// Return whether the ray is valid or not. + /// + /// It is possible for an invalid direction (zero length) + /// to be passed to the constructor. When this happens, + /// the constructor divides by zero, leaving Inf in all + /// components. + VTKM_EXEC_CONT + bool IsValid() const; + + /// Compute a point along the line. \a param values > 0 lie on the ray. + VTKM_EXEC_CONT + Vector Evaluate(CoordType param) const; + + /// Return the minmum distance from \a point to this line/ray. + /// + /// Note that when the direction has zero length, this simplifies + /// the distance between \a point and the ray's origin. + /// Otherwise, the distance returned is either the perpendicular + /// distance from \a point to the line or the distance + /// to the ray's origin. + VTKM_EXEC_CONT + CoordType DistanceTo(const Vector& point) const; + + /// Return the minimum distance between the ray/line and \a point. + VTKM_EXEC_CONT + CoordType DistanceTo(const Vector& point, CoordType& param, Vector& projectedPoint) const; + + /// Compute the non-degenerate point where two 2-D rays intersect, or return false. + /// + /// If true is returned, then the rays intersect in a unique point and + /// \a point is set to that location. + /// + /// If false is returned, then either + /// (1) the rays are parallel and may be + /// coincident (intersect everywhere) or offset (intersect nowhere); or + /// (2) the lines intersect but not the rays (because the the intersection + /// occurs in the negative parameter space of one or both rays). + /// In the latter case (2), the \a point is still set to the location of the + /// intersection. + /// + /// The tolerance \a tol is the minimum acceptable denominator used to + /// compute the intersection point coordinates and thus dictates the + /// maximum distance from the segments at which intersections will be + /// reported as valid. + template ::type = 0> + VTKM_EXEC_CONT bool Intersect(const Ray& other, + Vector& point, + CoordType tol = 0.f); +}; + +/// Represent a finite line segment with a pair of points +template +struct LineSegment +{ + static constexpr int Dimension = Dim; + using Vector = vtkm::Vec; + Vector Endpoints[2]; + + /// Construct a default segment from (0,0) to (1,0). + template ::type = 0> + VTKM_EXEC_CONT LineSegment(); + + /// Construct a default segment from (0,0,0) to (1,0,0). + template ::type = 0> + VTKM_EXEC_CONT LineSegment(); + + /// Construct a segment spanning points \a p0 and \a p1. + VTKM_EXEC_CONT + LineSegment(const Vector& p0, const Vector& p1); + + /// Return whether this line segment has an infinitesimal extent (i.e., whether the endpoints are coincident) + VTKM_EXEC_CONT + bool IsSingular(CoordType tol2 = static_cast(1.0e-6f)) const; + + /// Construct a plane bisecting this line segment (only when Dimension is 3). + template ::type = 0> + VTKM_EXEC_CONT Plane PerpendicularBisector() const; + + /// Construct a perpendicular bisector to this line segment (only when Dimension is 2). + template ::type = 0> + VTKM_EXEC_CONT Ray PerpendicularBisector() const; + + /// Return the midpoint of the line segment + VTKM_EXEC_CONT + Vector Center() const { return this->Evaluate(0.5f); } + + /// Return the vector pointing to endpoint 1 from endpoint 0. + /// This vector is not of unit length and in the case of + /// degenerate lines, may have zero length. + /// Call vtkm::Normal() on the return value if you want a normalized result. + VTKM_EXEC_CONT + Vector Direction() const { return this->Endpoints[1] - this->Endpoints[0]; } + + /// Compute a point along the line. \a param values in [0,1] lie on the line segment. + VTKM_EXEC_CONT + Vector Evaluate(CoordType param) const; + + /// Return the minmum distance from \a point to this line segment. + /// + /// Note that when the endpoints are coincident, this simplifies + /// the distance between \a point and either endpoint. + /// Otherwise, the distance returned is either the perpendicular + /// distance from \a point to the line segment or the distance + /// to the nearest endpoint (whichever is smaller). + VTKM_EXEC_CONT + CoordType DistanceTo(const Vector& point) const; + + /// Return the minimum distance between the line segment and \a point. + VTKM_EXEC_CONT + CoordType DistanceTo(const Vector& point, CoordType& param, Vector& projectedPoint) const; + + /// Compute the non-degenerate point where two (infinite) 2-D line segments intersect, or return false. + /// + /// If true is returned, then the lines intersect in a unique point and + /// \a point is set to that location. + /// + /// If false is returned, then the lines are parallel and either they are + /// coincident (intersect everywhere) or offset (intersect nowhere). + /// + /// The tolerance \a tol is the minimum acceptable denominator used to + /// compute the intersection point coordinates and thus dictates the + /// maximum distance from the segments at which intersections will be + /// reported as valid. + template ::type = 0> + VTKM_EXEC_CONT bool IntersectInfinite(const LineSegment& other, + Vector& point, + CoordType tol = 0.f); +}; + +/// Represent a plane with a base point (origin) and normal vector. +template +struct Plane +{ + using Vector = vtkm::Vec; + Vector Origin; + Vector Normal; + + /// Construct a default plane whose base point is the origin and whose normal is (0,0,1) + VTKM_EXEC_CONT + Plane(); + + /// Construct a plane with the given \a origin and \a normal. + VTKM_EXEC_CONT + Plane(const Vector& origin, const Vector& normal, CoordType tol2 = static_cast(1e-8f)); + + /// Return true if the plane's normal is well-defined to within the given tolerance. + VTKM_EXEC_CONT + bool IsValid() const { return !vtkm::IsInf(this->Normal[0]); } + + /// Return the **signed** distance from the plane to the point. + VTKM_EXEC_CONT + CoordType DistanceTo(const Vector& point) const; + + /// Return the closest point in the plane to the given point. + VTKM_EXEC_CONT + Vector ClosestPoint(const Vector& point) const; + + /// Intersect this plane with the ray (or line if the ray is two-sided). + /// + /// Returns true if there is a non-degenrate intersection (i.e., an isolated point of intersection). + /// Returns false if there is no intersection *or* if the intersection is degenerate (i.e., the + /// entire ray/line lies in the plane). + /// In the latter case, \a lineInPlane will be true upon exit. + /// + /// If this method returns true, then \a parameter will be set to a number indicating + /// where along the ray/line the plane hits and \a point will be set to that location. + /// If the input is a ray, the \a parameter will be non-negative. + template + VTKM_EXEC_CONT bool Intersect(const Ray& ray, + CoordType& parameter, + Vector& point, + bool& lineInPlane, + CoordType tol = CoordType(1e-6f)) const; + + /// Intersect this plane with the line \a segment. + /// + /// Returns true if there is a non-degenrate intersection (i.e., an isolated point of intersection). + /// Returns false if there is no intersection *or* if the intersection is degenerate (i.e., the + /// entire line segment lies in the plane). + /// In the latter case, \a lineInPlane will be true upon exit. + /// + /// If this method returns true, then \a parameter will be set to a number in [0,1] indicating + /// where along the line segment the plane hits. + VTKM_EXEC_CONT + bool Intersect(const LineSegment& segment, + CoordType& parameter, + bool& lineInPlane) const; + + /// Intersect this plane with the line \a segment. + /// + /// Returns true if there is a non-degenrate intersection (i.e., an isolated point of intersection). + /// Returns false if there is no intersection *or* if the intersection is degenerate (i.e., the + /// entire line segment lines in the plane). + /// In the latter case, \a lineInPlane will be true upon exit. + /// + /// If this method returns true, then \a parameter will be set to a number in [0,1] indicating + /// where along the line segment the plane hits and \a point will be set to that location. + VTKM_EXEC_CONT + bool Intersect(const LineSegment& segment, + CoordType& parameter, + Vector& point, + bool& lineInPlane) const; + + /// Intersect this plane with another plane. + /// + /// Returns true if there is a non-degenrate intersection (i.e., a line of intersection). + /// Returns false if there is no intersection *or* if the intersection is degenerate + /// (i.e., the planes are coincident). + /// In the latter case, \a coincident will be true upon exit and \a segment will + /// unmodified. + /// + /// If this method returns true, then the resulting \a segment will have its + /// base point on the line of intersection and its second point will be a unit + /// length away in the direction of the cross produce of the input plane normals + /// (this plane crossed with the \a other). + /// + /// The tolerance \a tol is the minimum squared length of the cross-product + /// of the two plane normals. It is also compared to the squared distance of + /// the base point of \a other away from \a this plane when considering whether + /// the planes are coincident. + VTKM_EXEC_CONT + bool Intersect(const Plane& other, + Ray& ray, + bool& coincident, + CoordType tol2 = static_cast(1e-6f)) const; +}; + +/// Represent a sphere of the given \a Dimension. +/// If a constructor is given an invalid specification, then +/// the Radius of the resulting sphere will be -1. +template +struct Sphere +{ + static constexpr int Dimension = Dim; + using Vector = vtkm::Vec; + Vector Center; + CoordType Radius; + + /// Construct a default sphere (unit radius at the origin). + VTKM_EXEC_CONT Sphere(); + + /// Construct a sphere from a center point and radius. + VTKM_EXEC_CONT Sphere(const Vector& center, CoordType radius); + + /// Return true if the sphere is valid (i.e., has a strictly positive radius). + VTKM_EXEC_CONT + bool IsValid() const { return this->Radius > 0.f; } + + /// Return whether the point lies strictly inside the sphere. + VTKM_EXEC_CONT + bool Contains(const Vector& point, CoordType tol2 = 0.f) const; + + /// Classify a point as inside (-1), on (0), or outside (+1) of the sphere. + /// + /// The tolerance \a tol2 is the maximum allowable difference in squared + /// magnitude between the squared radius and the squared distance between + /// the \a point and Center. + VTKM_EXEC_CONT + int Classify(const Vector& point, CoordType tol2 = 0.f) const; +}; + +// ----------------------------------------------------------------------------- +// Synonyms +// +// These "using" statements aim to make it easier to use the templated +// structs above when working with a particular dimension and/or the +// default floating-point type. + +/// Lines are two-sided rays: +template +using Line = Ray; + +// Shortcuts for 2D and 3D rays, lines, and line segments: +template +using Ray2 = Ray; +template +using Ray3 = Ray; +template +using Line2 = Line; +template +using Line3 = Line; +template +using LineSegment2 = LineSegment; +template +using LineSegment3 = LineSegment; + +/// Circle is an alias for a 2-Dimensional sphere. +template +using Circle = Sphere; + +// Aliases for d-dimensional spheres. +template +using Sphere2 = Sphere; +template +using Sphere3 = Sphere; + +// Shortcuts for default floating-point types +using Ray2d = Ray2; +using Ray3d = Ray3; +using Line2d = Line2; +using Line3d = Line3; +using LineSegment2d = LineSegment2; +using LineSegment3d = LineSegment3; +using Plane3d = Plane; +using Circle2d = Circle; +using Sphere2d = Sphere2; +using Sphere3d = Sphere3; + +// ----------------------------------------------------------------------------- +// Construction techniques +// +// These are free functions that create instances of geometric structs by taking +// in data that is not identical to the state of the struct and converting it +// into state for the struct. + +/// Construct a plane from a point plus one of: a line, a ray, or a line segment. +/// +/// The plane returned will contain the point and the line/ray/segment. +/// The plane normal will be the cross product of the line/ray/segment's direction +/// and the vector from the line/ray/segment's origin to the given \a point. +/// If the \a point is collinear with the line/ray/line-segment, an invalid +/// plane will be returned. +template +VTKM_EXEC_CONT vtkm::Plane make_PlaneFromPointAndLine( + const vtkm::Vec& point, + const vtkm::Ray& ray, + CoordType tol2 = static_cast(1e-8f)); + +template +VTKM_EXEC_CONT vtkm::Plane make_PlaneFromPointAndLineSegment( + const vtkm::Vec& point, + const vtkm::LineSegment3& segment, + CoordType tol2 = static_cast(1e-8f)); + +/// Construct a circle from 3 points. +template +VTKM_EXEC_CONT vtkm::Circle make_CircleFrom3Points( + const typename vtkm::Vec& p0, + const typename vtkm::Vec& p1, + const typename vtkm::Vec& p2, + CoordType tol = static_cast(1e-6f)); + +/// Construct a sphere from 4 points. +template +VTKM_EXEC_CONT vtkm::Sphere make_SphereFrom4Points( + const vtkm::Vec& a0, + const vtkm::Vec& a1, + const vtkm::Vec& a2, + const vtkm::Vec& a3, + CoordType tol = static_cast(1e-6f)); + +} // namespace vtkm + +#include "vtkm/Geometry.hxx" + +#endif // vtk_m_Geometry_h diff --git a/vtkm/Geometry.hxx b/vtkm/Geometry.hxx new file mode 100644 index 000000000..cb91818e7 --- /dev/null +++ b/vtkm/Geometry.hxx @@ -0,0 +1,647 @@ +//============================================================================ +// 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 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +// Copyright 2014 UT-Battelle, LLC. +// Copyright 2014 Los Alamos National Security. +// +// Under the terms of Contract DE-NA0003525 with NTESS, +// 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. +//============================================================================ + +namespace vtkm +{ + +// ----------------------------------------------------------------------------- +// Ray + +template +template ::type> +Ray::Ray() + : Origin{ 0.f } + , Direction{ 1.f, 0.f } +{ +} + +template +template ::type> +Ray::Ray() + : Origin{ 0.f } + , Direction{ 1.f, 0.f, 0.f } +{ +} + +template +Ray::Ray(const LineSegment& segment) + : Origin(segment.Endpoints[0]) + , Direction(vtkm::Normal(segment.Direction())) +{ +} + +template +Ray::Ray(const Vector& point, const Vector& direction) + : Origin(point) + , Direction(vtkm::Normal(direction)) +{ +} + +template +typename Ray::Vector Ray::Evaluate( + CoordType param) const +{ + auto pointOnLine = this->Origin + this->Direction * param; + return pointOnLine; +} + +template +bool Ray::IsValid() const +{ +// At least on Ubuntu 17.10, cuda 9.1 will fail with an internal +// compiler error when calling vtkm::IsInf() here. But the fix +// below works. The fix should be removed as soon as our dashboards +// allow it. +#if __CUDACC_VER_MAJOR__ == 9 && __CUDACC_VER_MINOR__ == 1 + return !isinf(this->Direction[0]); +#else + return !vtkm::IsInf(this->Direction[0]); +#endif +} + +template +CoordType Ray::DistanceTo(const Vector& point) const +{ + Vector closest; + CoordType param; + return this->DistanceTo(point, param, closest); +} + +template +CoordType Ray::DistanceTo(const Vector& point, + CoordType& param, + Vector& projectedPoint) const +{ + const auto& dir = this->Direction; + auto mag2 = vtkm::MagnitudeSquared(dir); + if (mag2 <= static_cast(0.0f)) + { + // We have a point, not a line segment. + projectedPoint = this->Origin; + param = static_cast(0.0f); + return vtkm::Magnitude(point - this->Origin); + } + + // Find the closest point on the line, then clamp to the ray if the + // parameter value is negative. + param = vtkm::Dot(point - this->Origin, dir) / mag2; + if (!TwoSided) + { + param = vtkm::Max(param, static_cast(0.0f)); + } + + // Compute the distance between the closest point and the input point. + projectedPoint = this->Evaluate(param); + auto dist = vtkm::Magnitude(point - projectedPoint); + return dist; +} + +template +template ::type> +bool Ray::Intersect(const Ray& other, + Vector& point, + CoordType tol) +{ + auto d1 = this->Direction; + auto d2 = other.Direction; + auto denom = d1[0] * d2[1] - d1[1] * d2[0]; + if (vtkm::Abs(denom) < tol) + { // The lines are coincident or at least parallel. + return false; + } + const auto& a = this->Origin; + const auto& b = other.Origin; + CoordType numerU = a[1] * d2[0] + d2[1] * b[0] - b[1] * d2[0] - d2[1] * a[0]; + + CoordType uParam = numerU / denom; + point = a + uParam * d1; + if (TwoSided && OtherTwoSided) + { + return true; + } + else + { + CoordType numerV = d1[0] * (a[1] - b[1]) - d1[1] * (a[0] - b[0]); + CoordType vParam = numerV / denom; + return (TwoSided || (uParam + tol) > 0) && (OtherTwoSided || (vParam + tol) > 0); + } +} + +// ----------------------------------------------------------------------------- +// LineSegment + +template +template ::type> +LineSegment::LineSegment() + : Endpoints{ { 0.f }, { 1.f, 0.f } } +{ +} + +template +template ::type> +LineSegment::LineSegment() + : Endpoints{ { 0.f }, { 1.f, 0.f, 0.f } } +{ +} + +template +LineSegment::LineSegment(const Vector& p0, const Vector& p1) + : Endpoints{ p0, p1 } +{ +} + +template +bool LineSegment::IsSingular(CoordType tol2) const +{ + return vtkm::MagnitudeSquared(this->Direction()) < tol2; +} + +template +template ::type> +Ray LineSegment::PerpendicularBisector() const +{ + const Vector dir = this->Direction(); + const Vector perp(-dir[1], dir[0]); + const Vector mid = this->Center(); + return Ray(mid, perp); +} + +template +template ::type> +Plane LineSegment::PerpendicularBisector() const +{ + return Plane(this->Center(), this->Direction()); +} + +template +typename LineSegment::Vector LineSegment::Evaluate( + CoordType param) const +{ + auto pointOnLine = this->Endpoints[0] * (1.0f - param) + this->Endpoints[1] * param; + return pointOnLine; +} + +template +CoordType LineSegment::DistanceTo(const Vector& point) const +{ + Vector closest; + CoordType param; + return this->DistanceTo(point, param, closest); +} + +template +CoordType LineSegment::DistanceTo(const Vector& point, + CoordType& param, + Vector& projectedPoint) const +{ + auto dir = this->Endpoints[1] - this->Endpoints[0]; + auto mag2 = vtkm::MagnitudeSquared(dir); + if (mag2 <= static_cast(0.0f)) + { + // We have a point, not a line segment. + projectedPoint = this->Endpoints[0]; + param = static_cast(0.0f); + return vtkm::Magnitude(point - this->Endpoints[0]); + } + + // Find the closest point on the line, then clamp to the line segment + param = vtkm::Clamp(vtkm::Dot(point - this->Endpoints[0], dir) / mag2, + static_cast(0.0f), + static_cast(1.0f)); + + // Compute the distance between the closest point and the input point. + projectedPoint = this->Evaluate(param); + auto dist = vtkm::Magnitude(point - projectedPoint); + return dist; +} + +template +template ::type> +bool LineSegment::IntersectInfinite(const LineSegment& other, + Vector& point, + CoordType tol) +{ + auto d1 = this->Direction(); + auto d2 = other.Direction(); + auto denom = d1[0] * d2[1] - d1[1] * d2[0]; + if (vtkm::Abs(denom) < tol) + { // The lines are coincident or at least parallel. + return false; + } + const auto& a = this->Endpoints; + const auto& b = other.Endpoints; + CoordType numerX = (a[0][0] * a[1][1] - a[0][1] * a[1][0]) * -d2[0] - + (b[0][0] * b[1][1] - b[0][1] * b[1][0]) * -d1[0]; + CoordType numerY = (a[0][0] * a[1][1] - a[0][1] * a[1][0]) * -d2[1] - + (b[0][0] * b[1][1] - b[0][1] * b[1][0]) * -d1[1]; + point = Vector(numerX / denom, numerY / denom); + return true; +} + +// ----------------------------------------------------------------------------- +// Plane + +template +Plane::Plane() + : Origin{ 0.f, 0.f, 0.f } + , Normal{ 0.f, 0.f, 1.f } +{ +} + +template +Plane::Plane(const Vector& origin, const Vector& normal, CoordType tol2) + : Origin(origin) + , Normal(vtkm::Normal(normal)) +{ + if (tol2 > 0.0f && vtkm::MagnitudeSquared(normal) < tol2) + { + auto inf = vtkm::Infinity(); + Normal = vtkm::Vec(inf, inf, inf); + } +} + +template +CoordType Plane::DistanceTo(const Vector& point) const +{ + auto dist = vtkm::Dot(point - this->Origin, this->Normal); + return dist; +} + +template +typename Plane::Vector Plane::ClosestPoint(const Vector& point) const +{ + auto vop = vtkm::Project(point - this->Origin, this->Normal); + auto closest = point - vop; + return closest; +} + +template +template +bool Plane::Intersect(const Ray& ray, + CoordType& parameter, + Vector& point, + bool& lineInPlane, + CoordType tol) const +{ + CoordType d0 = this->DistanceTo(ray.Origin); + CoordType dirDot = vtkm::Dot(this->Normal, ray.Direction); + // If the ray/line lies parallel to the plane, the intersection is degenerate: + if (vtkm::Abs(dirDot) < tol) + { + lineInPlane = (vtkm::Abs(d0) < tol); + return false; + } + lineInPlane = false; + parameter = -d0 / dirDot; + // If we have a ray (not a line) and it points away from the + // side of the plane where its origin lies, then there is no + // intersection. + if (!IsTwoSided) + { + if (parameter < 0.0f) + { + return false; + } + } + + // Check whether an endpoint lies in the plane: + if (vtkm::Abs(d0) < tol) + { + parameter = static_cast(0.0f); + point = ray.Origin; + return true; + } + + // The perpendicular distance of the origin to the plane + // forms one side of a triangle whose hypotenuse is the + // parameter value (because ray.Direction has unit length). + // The dot product of the plane normal and ray direction + // is the cosine of the angle between the hypotenuse and + // the shortest path to the plane, so.... + point = ray.Origin + parameter * ray.Direction; + return true; +} + +template +bool Plane::Intersect(const LineSegment& segment, + CoordType& parameter, + bool& lineInPlane) const +{ + Vector point; + return this->Intersect(segment, parameter, point, lineInPlane); +} + +template +bool Plane::Intersect(const LineSegment& segment, + CoordType& parameter, + Vector& point, + bool& lineInPlane) const +{ + CoordType d0 = this->DistanceTo(segment.Endpoints[0]); + CoordType d1 = this->DistanceTo(segment.Endpoints[1]); + if (d0 == 0 && d1 == 0) + { + // The entire segment lies in the plane. + lineInPlane = true; + return true; + } + + lineInPlane = false; + // Check whether an endpoint lies in the plane: + if (d0 == 0) + { + parameter = static_cast(0.0f); + point = segment.Endpoints[0]; + return true; + } + if (d1 == 0) + { + parameter = static_cast(1.0f); + point = segment.Endpoints[1]; + return true; + } + + // See whether endpoints lie on opposite sides of the plane: + // + // Note that the vtkm::SignBit comments below cause an internal compiler + // error on cuda 9.1, ubuntu 17.10 or they would be used: + bool c0 = d0 < 0; // !!vtkm::SignBit(d0); + bool c1 = d1 < 0; // !!vtkm::SignBit(d1); + CoordType a0 = vtkm::Abs(d0); + CoordType a1 = vtkm::Abs(d1); + if (c0 == c1) + { + // Both endpoints lie to the same side of the plane, so + // there is no intersection. Report the closest endpoint. + parameter = static_cast(a0 < a1 ? 0.0f : 1.0f); + point = segment.Endpoints[a0 < a1 ? 0 : 1]; + return false; + } + + // Endpoint distances have the opposite sign; there must be + // an intersection. It must occur at distance 0, and distance + // varies linearly from d0 to d1, so... + parameter = a0 / (a0 + a1); + point = segment.Endpoints[0] * (1.0f - parameter) + segment.Endpoints[1] * parameter; + return true; +} + +template +bool Plane::Intersect(const Plane& other, + Ray& ray, + bool& coincident, + CoordType tol2) const +{ + auto dir = vtkm::Cross(this->Normal, other.Normal); + auto mag2 = vtkm::MagnitudeSquared(dir); + if (mag2 < tol2) + { // The planes are parallel. + auto dist = this->DistanceTo(other.Origin); + coincident = dist * dist < tol2; + return false; + } + // The planes intersect. We want to find a point on the new plane + // and we want it to be near the other plane base points to avoid + // precision issues in the future. + // So, project each plane origin to the other plane along a line + // perpendicular to the plane and the output line. Both of these + // points are on the output line. Average the two points. The result + // will still be on the line and will be closer to the two base points. + auto nn = vtkm::Normal(dir); + auto moveDir01 = vtkm::Cross(this->Normal, nn); + auto moveDir02 = vtkm::Cross(other.Normal, nn); + Ray bra(this->Origin, moveDir01); + Ray brb(other.Origin, moveDir02); + vtkm::Vec p0a; + vtkm::Vec p0b; + CoordType pDummyA, pDummyB; + bool bDummyA, bDummyB; + auto tol = vtkm::Sqrt(tol2); + this->Intersect(brb, pDummyA, p0a, bDummyA, tol); + other.Intersect(bra, pDummyB, p0b, bDummyB, tol); + ray = vtkm::Ray((p0a + p0b) * 0.5f, nn); + return true; +} + +// ----------------------------------------------------------------------------- +// Sphere + +template +Sphere::Sphere() + : Center{ 0.f } + , Radius(static_cast(1.f)) +{ +} + +template +Sphere::Sphere(const Vector& center, CoordType radius) + : Center(center) + , Radius(radius <= 0.f ? static_cast(-1.0f) : radius) +{ +} + +template +bool Sphere::Contains(const Vector& point, CoordType tol2) const +{ + return this->Classify(point, tol2) < 0; +} + +template +int Sphere::Classify(const Vector& point, CoordType tol2) const +{ + if (!this->IsValid()) + { + return 1; // All points are outside invalid spheres. + } + auto d2 = vtkm::MagnitudeSquared(point - this->Center); + auto r2 = this->Radius * this->Radius; + return d2 < r2 - tol2 ? -1 : (d2 > r2 + tol2 ? +1 : 0); +} + +// ----------------------------------------------------------------------------- +// Construction techniques + +template +vtkm::Plane make_PlaneFromPointAndLine(const vtkm::Vec& point, + const vtkm::Ray& ray, + CoordType tol2) +{ + auto tmpDir = point - ray.Origin; + return vtkm::Plane(point, vtkm::Cross(ray.Direction, tmpDir), tol2); +} + +template +vtkm::Plane make_PlaneFromPointAndLineSegment( + const vtkm::Vec& point, + const vtkm::LineSegment3& segment, + CoordType tol2) +{ + auto tmpDir = point - segment.Endpoints[0]; + return vtkm::Plane(point, vtkm::Cross(segment.Direction(), tmpDir), tol2); +} + +template +vtkm::Circle make_CircleFrom3Points(const typename vtkm::Vec& p0, + const typename vtkm::Vec& p1, + const typename vtkm::Vec& p2, + CoordType tol) +{ + constexpr int Dim = 2; + using Vector = typename vtkm::Circle::Vector; + + vtkm::LineSegment l01(p0, p1); + vtkm::LineSegment l02(p0, p2); + auto pb01 = l01.PerpendicularBisector(); + auto pb02 = l02.PerpendicularBisector(); + Vector center; + CoordType radius; + if (!pb01.IsValid() || !pb02.IsValid()) + { + center = Vector(0.f, 0.f); + radius = -1.f; + return vtkm::Circle(center, radius); + } + auto isValid = pb01.Intersect(pb02, center, tol); + radius = isValid ? vtkm::Magnitude(center - p0) : static_cast(-1.0f); + if (!isValid) + { // Initialize center to something. + center = Vector(vtkm::Nan(), vtkm::Nan()); + } + return vtkm::Circle(center, radius); +} + +template +vtkm::Sphere make_SphereFrom4Points(const vtkm::Vec& a0, + const vtkm::Vec& a1, + const vtkm::Vec& a2, + const vtkm::Vec& a3, + CoordType tol) +{ + // Choose p3 such that the min(p3 - p[012]) is larger than any other choice of p3. + // From: http://steve.hollasch.net/cgindex/geometry/sphere4pts.html, + // retrieved 2018-06-12 and paraphrased in terms of local variable names: + // + // If circlePointInPlaneOfP3-p3 is much smaller than + // circlePointInPlaneOfP3-circleCenterWorld, then the + // sphere center will be very close to circleCenterWorld + // and subject to error. + // It's best to choose p3 so that the least of p0-p3, p1-p3, and + // p2-p3 is larger than for any other. + + CoordType d0 = vtkm::MagnitudeSquared(a1 - a0); + CoordType d1 = vtkm::MagnitudeSquared(a2 - a0); + CoordType d2 = vtkm::MagnitudeSquared(a3 - a0); + CoordType d3 = vtkm::MagnitudeSquared(a2 - a1); + CoordType d4 = vtkm::MagnitudeSquared(a3 - a1); + CoordType d5 = vtkm::MagnitudeSquared(a3 - a2); + CoordType sel0 = vtkm::Min(d0, vtkm::Min(d1, d2)); + CoordType sel1 = vtkm::Min(d0, vtkm::Min(d3, d4)); + CoordType sel2 = vtkm::Min(d1, vtkm::Min(d3, d5)); + CoordType sel3 = vtkm::Min(d2, vtkm::Min(d4, d5)); + CoordType selm = vtkm::Max(vtkm::Max(sel0, sel1), vtkm::Max(sel2, sel3)); + + vtkm::Vec p0 = a0; + vtkm::Vec p1 = a1; + vtkm::Vec p2 = a2; + vtkm::Vec p3 = a3; + if (sel0 == selm) + { + p3 = a0; + p0 = a3; + } + else if (sel1 == selm) + { + p3 = a1; + p1 = a3; + } + else if (sel2 == selm) + { + p3 = a2; + p2 = a3; + } + else // sel3 == selm + { + // do nothing. + } + + vtkm::Vec, 3> axes; + axes[1] = p1 - p0; + axes[2] = p2 - p0; + axes[0] = vtkm::Cross(axes[1], axes[2]); + vtkm::Vec, 3> basis; + int rank = vtkm::Orthonormalize(axes, basis, tol); + if (rank < 3) + { + vtkm::Sphere invalid; + invalid.Radius = -1.0f; + return invalid; + } + + // Project points to plane and fit a circle. + auto p0_p = vtkm::Vec{ 0.f }; // This is p0's new coordinate... + auto p1_p = vtkm::Vec(vtkm::ProjectedDistance(axes[1], basis[1]), + vtkm::ProjectedDistance(axes[1], basis[2])); + auto p2_p = vtkm::Vec(vtkm::ProjectedDistance(axes[2], basis[1]), + vtkm::ProjectedDistance(axes[2], basis[2])); + + auto circle = make_CircleFrom3Points(p0_p, p1_p, p2_p); + if (!circle.IsValid()) + { + vtkm::Sphere invalid; + invalid.Radius = -1.0f; + return invalid; + } + + vtkm::Vec circleCenterWorld = + p0 + circle.Center[0] * basis[1] + circle.Center[1] * basis[2]; + + vtkm::Line3 centerRay(circleCenterWorld, basis[0]); + // If our remaining unused point (p3) lines on centerRay, + // use one of the other points to locate the sphere's center: + vtkm::Vec circlePointInPlaneOfP3; + if (vtkm::Abs(centerRay.DistanceTo(p3)) < tol) + { + circlePointInPlaneOfP3 = p0; + } + else + { + Plane pp3(circleCenterWorld, basis[0]); + circlePointInPlaneOfP3 = + circleCenterWorld + vtkm::Normal(pp3.ClosestPoint(p3) - circleCenterWorld) * circle.Radius; + } + + auto bisectorPlane = + vtkm::LineSegment3(circlePointInPlaneOfP3, p3).PerpendicularBisector(); + vtkm::Vec sphereCenter; + CoordType param; + bool lineInPlane; + if (!bisectorPlane.Intersect(centerRay, param, sphereCenter, lineInPlane, tol)) + { + vtkm::Sphere invalid; + invalid.Radius = -1.0f; + return invalid; + } + auto sphereRadius = vtkm::Magnitude(sphereCenter - p3); + vtkm::Sphere result(sphereCenter, sphereRadius); + ; + return result; +} + +} // namespace vtkm diff --git a/vtkm/Math.h b/vtkm/Math.h index 224981db2..acb74a951 100644 --- a/vtkm/Math.h +++ b/vtkm/Math.h @@ -1815,6 +1815,18 @@ static inline VTKM_EXEC_CONT T Min(const T& x, const T& y) return detail::Min(x, y, typename vtkm::TypeTraits::DimensionalityTag()); } +/// Clamp \p x to the given range. +/// + +inline VTKM_EXEC_CONT vtkm::Float32 Clamp(vtkm::Float32 x, vtkm::Float32 lo, vtkm::Float32 hi) +{ + return x > lo ? (x < hi ? x : hi) : lo; +} + +inline VTKM_EXEC_CONT vtkm::Float64 Clamp(vtkm::Float64 x, vtkm::Float64 lo, vtkm::Float64 hi) +{ + return x > lo ? (x < hi ? x : hi) : lo; +} //----------------------------------------------------------------------------- //#ifdef VTKM_CUDA diff --git a/vtkm/Math.h.in b/vtkm/Math.h.in index 26492cdc1..af810f8d3 100644 --- a/vtkm/Math.h.in +++ b/vtkm/Math.h.in @@ -637,6 +637,18 @@ static inline VTKM_EXEC_CONT T Min(const T& x, const T& y) return detail::Min(x, y, typename vtkm::TypeTraits::DimensionalityTag()); } +/// Clamp \p x to the given range. +/// + +inline VTKM_EXEC_CONT vtkm::Float32 Clamp(vtkm::Float32 x, vtkm::Float32 lo, vtkm::Float32 hi) +{ + return x > lo ? (x < hi ? x : hi) : lo; +} + +inline VTKM_EXEC_CONT vtkm::Float64 Clamp(vtkm::Float64 x, vtkm::Float64 lo, vtkm::Float64 hi) +{ + return x > lo ? (x < hi ? x : hi) : lo; +} //----------------------------------------------------------------------------- //#ifdef VTKM_CUDA diff --git a/vtkm/TypeListTag.h b/vtkm/TypeListTag.h index b49e3157e..057ade54e 100644 --- a/vtkm/TypeListTag.h +++ b/vtkm/TypeListTag.h @@ -92,6 +92,19 @@ struct TypeListTagFieldVec4 { }; +/// A list containing common types for floating-point vectors. Specifically contains +/// floating point vectors of size 2, 3, and 4 with floating point components. +/// Scalars are not included. +/// +struct TypeListTagFloatVec : vtkm::ListTagBase, + vtkm::Vec, + vtkm::Vec, + vtkm::Vec, + vtkm::Vec, + vtkm::Vec> +{ +}; + /// A list containing common types for values in fields. Specifically contains /// floating point scalars and vectors of size 2, 3, and 4 with floating point /// components. diff --git a/vtkm/VectorAnalysis.h b/vtkm/VectorAnalysis.h index 939b57d27..978e5c4d2 100644 --- a/vtkm/VectorAnalysis.h +++ b/vtkm/VectorAnalysis.h @@ -208,6 +208,84 @@ TriangleNormal(const vtkm::Vec& a, const vtkm::Vec& b, const vtkm::V return vtkm::Cross(b - a, c - a); } +//----------------------------------------------------------------------------- +/// \brief Project a vector onto another vector. +/// +/// This method computes the orthogonal projection of the vector v onto u; +/// that is, it projects its first argument onto its second. +/// +/// Note that if the vector \a u has zero length, the output +/// vector will have all its entries equal to NaN. +template +VTKM_EXEC_CONT vtkm::Vec Project(const vtkm::Vec& v, const vtkm::Vec& u) +{ + T uu = vtkm::Dot(u, u); + T uv = vtkm::Dot(u, v); + T factor = uv / uu; + vtkm::Vec result = factor * u; + return result; +} + +//----------------------------------------------------------------------------- +/// \brief Project a vector onto another vector, returning only the projected distance. +/// +/// This method computes the orthogonal projection of the vector v onto u; +/// that is, it projects its first argument onto its second. +/// +/// Note that if the vector \a u has zero length, the output will be NaN. +template +VTKM_EXEC_CONT T ProjectedDistance(const vtkm::Vec& v, const vtkm::Vec& u) +{ + T uu = vtkm::Dot(u, u); + T uv = vtkm::Dot(u, v); + T factor = uv / uu; + return factor; +} + +//----------------------------------------------------------------------------- +/// \brief Perform Gram-Schmidt orthonormalization for 3-D vectors. +/// +/// See https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process for details. +/// The first output vector will always be parallel to the first input vector. +/// The remaining output vectors will be orthogonal and unit length and have +/// the same handedness as their corresponding input vectors. +/// +/// This method is geometric. +/// It does not require a matrix solver. +/// However, unlike the algebraic eigensolver techniques which do use matrix +/// inversion, this method may return zero-length output vectors if some input +/// vectors are collinear. The number of non-zero (to within the specified +/// tolerance, \a tol ) output vectors is the return value. +template +VTKM_EXEC_CONT int Orthonormalize(const vtkm::Vec, N>& inputs, + vtkm::Vec, N>& outputs, + T tol = static_cast(1e-6)) +{ + int j = 0; // j is the number of non-zero-length, non-collinear inputs encountered. + vtkm::Vec, N> u; + for (int i = 0; i < N; ++i) + { + u[j] = inputs[i]; + for (int k = 0; k < j; ++k) + { + u[j] -= vtkm::Project(inputs[i], u[k]); + } + T rmag = vtkm::RMagnitude(u[j]); + if (rmag * tol > 1.0) + { + // skip this vector, it is zero-length or collinear with others. + continue; + } + outputs[j] = rmag * u[j]; + ++j; + } + for (int i = j; i < N; ++i) + { + outputs[j] = Vec{ 0. }; + } + return j; +} + } // namespace vtkm #endif //vtk_m_VectorAnalysis_h diff --git a/vtkm/cont/cuda/testing/CMakeLists.txt b/vtkm/cont/cuda/testing/CMakeLists.txt index 8b6473e1a..d5a4e01f4 100644 --- a/vtkm/cont/cuda/testing/CMakeLists.txt +++ b/vtkm/cont/cuda/testing/CMakeLists.txt @@ -28,6 +28,7 @@ set(unit_tests UnitTestCudaDataSetExplicit.cu UnitTestCudaDataSetSingleType.cu UnitTestCudaDeviceAdapter.cu + UnitTestCudaGeometry.cu UnitTestCudaImplicitFunction.cu UnitTestCudaMath.cu UnitTestCudaShareUserProvidedManagedMemory.cu diff --git a/vtkm/cont/cuda/testing/UnitTestCudaGeometry.cu b/vtkm/cont/cuda/testing/UnitTestCudaGeometry.cu new file mode 100644 index 000000000..d5caa0782 --- /dev/null +++ b/vtkm/cont/cuda/testing/UnitTestCudaGeometry.cu @@ -0,0 +1,38 @@ +//============================================================================= +// +// 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 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +// Copyright 2015 UT-Battelle, LLC. +// Copyright 2015 Los Alamos National Security. +// +// Under the terms of Contract DE-NA0003525 with NTESS, +// 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. +// +//============================================================================= + +// Make sure that the tested code is using the device adapter specified. This +// is important in the long run so we don't, for example, use the CUDA device +// for a part of an operation where the TBB device was specified. +#define VTKM_DEVICE_ADAPTER VTKM_DEVICE_ADAPTER_ERROR + +#include +#include +#include + +int UnitTestCudaGeometry(int, char* []) +{ + auto tracker = vtkm::cont::GetGlobalRuntimeDeviceTracker(); + tracker.ForceDevice(vtkm::cont::DeviceAdapterTagCuda{}); + return vtkm::cont::testing::Testing::Run( + UnitTestGeometryNamespace::RunGeometryTests); +} diff --git a/vtkm/cont/serial/testing/CMakeLists.txt b/vtkm/cont/serial/testing/CMakeLists.txt index ada4701c7..b1768dfd5 100644 --- a/vtkm/cont/serial/testing/CMakeLists.txt +++ b/vtkm/cont/serial/testing/CMakeLists.txt @@ -28,6 +28,7 @@ set(unit_tests UnitTestSerialDataSetExplicit.cxx UnitTestSerialDataSetSingleType.cxx UnitTestSerialDeviceAdapter.cxx + UnitTestSerialGeometry.cxx UnitTestSerialImplicitFunction.cxx UnitTestSerialPointLocatorUniformGrid.cxx UnitTestSerialVirtualObjectHandle.cxx diff --git a/vtkm/cont/serial/testing/UnitTestSerialGeometry.cxx b/vtkm/cont/serial/testing/UnitTestSerialGeometry.cxx new file mode 100644 index 000000000..7295dea6c --- /dev/null +++ b/vtkm/cont/serial/testing/UnitTestSerialGeometry.cxx @@ -0,0 +1,38 @@ +//============================================================================= +// +// 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 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +// Copyright 2015 UT-Battelle, LLC. +// Copyright 2015 Los Alamos National Security. +// +// Under the terms of Contract DE-NA0003525 with NTESS, +// 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. +// +//============================================================================= + +// Make sure that the tested code is using the device adapter specified. This +// is important in the long run so we don't, for example, use the CUDA device +// for a part of an operation where the TBB device was specified. +#define VTKM_DEVICE_ADAPTER VTKM_DEVICE_ADAPTER_ERROR + +#include +#include +#include + +int UnitTestSerialGeometry(int, char* []) +{ + auto tracker = vtkm::cont::GetGlobalRuntimeDeviceTracker(); + tracker.ForceDevice(vtkm::cont::DeviceAdapterTagSerial{}); + return vtkm::cont::testing::Testing::Run( + UnitTestGeometryNamespace::RunGeometryTests); +} diff --git a/vtkm/testing/CMakeLists.txt b/vtkm/testing/CMakeLists.txt index 8fead85c6..cc0b521bc 100644 --- a/vtkm/testing/CMakeLists.txt +++ b/vtkm/testing/CMakeLists.txt @@ -21,6 +21,7 @@ set(headers Testing.h TestingMath.h + TestingGeometry.h OptionParser.h VecTraitsTests.h ) diff --git a/vtkm/testing/TestingGeometry.h b/vtkm/testing/TestingGeometry.h new file mode 100644 index 000000000..11cc955d4 --- /dev/null +++ b/vtkm/testing/TestingGeometry.h @@ -0,0 +1,498 @@ +//============================================================================= +// +// 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 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +// Copyright 2015 UT-Battelle, LLC. +// Copyright 2015 Los Alamos National Security. +// +// Under the terms of Contract DE-NA0003525 with NTESS, +// 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_testing_TestingGeometry_h +#define vtk_m_testing_TestingGeometry_h + +#include +#include + +#include +#include + +#include + +#include + +#include + +#define VTKM_MATH_ASSERT(condition, message) \ + if (!(condition)) \ + { \ + this->RaiseError(message); \ + } + +//----------------------------------------------------------------------------- +namespace UnitTestGeometryNamespace +{ + +class Coords +{ +public: + static constexpr vtkm::IdComponent NUM_COORDS = 5; + + template + VTKM_EXEC_CONT vtkm::Vec EndpointList(vtkm::Int32 i) const + { + vtkm::Float64 coords[NUM_COORDS][3] = { + { 1.0, 0.0, 0.0 }, { 0.0, 1.0, 0.0 }, { -1.0, 0.0, 0.0 }, + { -2.0, 0.0, 0.0 }, { 0.0, -2.0, 0.0 }, + }; + vtkm::Int32 j = i % NUM_COORDS; + return vtkm::Vec( + static_cast(coords[j][0]), static_cast(coords[j][1]), static_cast(coords[j][2])); + } + + template + VTKM_EXEC_CONT vtkm::Vec ClosestToOriginList(vtkm::Int32 i) const + { + vtkm::Float64 coords[NUM_COORDS][3] = { + { 0.5, 0.5, 0.0 }, { -0.5, 0.5, 0.0 }, { -1.0, 0.0, 0.0 }, + { -1.0, -1.0, 0.0 }, { 0.8, -0.4, 0.0 }, + }; + vtkm::Int32 j = i % NUM_COORDS; + return vtkm::Vec( + static_cast(coords[j][0]), static_cast(coords[j][1]), static_cast(coords[j][2])); + } + + template + VTKM_EXEC_CONT T DistanceToOriginList(vtkm::Int32 i) const + { + vtkm::Float64 coords[NUM_COORDS] = { + 0.707107, 0.707107, 1.0, 1.41421, 0.894427, + }; + vtkm::Int32 j = i % NUM_COORDS; + return static_cast(coords[j]); + } +}; + +//----------------------------------------------------------------------------- + +template +struct RayTests : public vtkm::exec::FunctorBase +{ + VTKM_EXEC + void operator()(vtkm::Id) const + { + { + using V2 = vtkm::Vec; + using Ray2 = vtkm::Ray2; + + Ray2 ray0; + VTKM_MATH_ASSERT(test_equal(ray0.Origin, V2(0., 0.)), "Bad origin for default 2D ray ctor."); + VTKM_MATH_ASSERT(test_equal(ray0.Direction, V2(1., 0.)), + "Bad direction for default 2D ray ctor."); + + // Test intersection + Ray2 ray1(V2(-1., 0.), V2(+1., +1.)); + Ray2 ray2(V2(+1., 0.), V2(-1., +1.)); + V2 point; + bool didIntersect = ray1.Intersect(ray2, point); + VTKM_MATH_ASSERT(test_equal(didIntersect, true), "Ray-pair 1 should intersect."); + VTKM_MATH_ASSERT(test_equal(point, V2(0., 1.)), "Ray-pair 1 should intersect at (0,1)."); + + // Test non intersection + Ray2 ray3(V2(-1., 0.), V2(-1., -1.)); + Ray2 ray4(V2(+1., 0.), V2(+1., -1.)); + + didIntersect = ray1.Intersect(ray4, point); + VTKM_MATH_ASSERT(test_equal(didIntersect, false), "Ray-pair 2 should not intersect."); + VTKM_MATH_ASSERT(test_equal(point, V2(0., 1.)), "Ray-pair 2 should intersect at (0,1)."); + + didIntersect = ray3.Intersect(ray2, point); + VTKM_MATH_ASSERT(test_equal(didIntersect, false), "Ray-pair 3 should not intersect."); + VTKM_MATH_ASSERT(test_equal(point, V2(0., 1.)), "Ray-pair 3 should intersect at (0,1)."); + + didIntersect = ray3.Intersect(ray4, point); + VTKM_MATH_ASSERT(test_equal(didIntersect, false), "Ray-pair 4 should not intersect."); + VTKM_MATH_ASSERT(test_equal(point, V2(0., 1.)), "Ray-pair 4 should intersect at (0,1)."); + } + + { + using V3 = vtkm::Vec; + + vtkm::Ray ray0; + VTKM_MATH_ASSERT(test_equal(ray0.Origin, V3(0., 0., 0.)), + "Bad origin for default 3D ray ctor."); + VTKM_MATH_ASSERT(test_equal(ray0.Direction, V3(1., 0., 0.)), + "Bad direction for default 3D ray ctor."); + } + } +}; + +template +struct TryRayTests +{ + template + void operator()(const T&) const + { + vtkm::cont::DeviceAdapterAlgorithm::Schedule(RayTests(), 1); + } +}; + +//----------------------------------------------------------------------------- + +template +struct LineSegmentTests : public vtkm::exec::FunctorBase +{ + VTKM_EXEC + void operator()(vtkm::Id) const + { + { + using V2 = vtkm::Vec; + using Line2 = vtkm::Line2; + + vtkm::LineSegment seg0; + VTKM_MATH_ASSERT(test_equal(seg0.Endpoints[0], V2(0., 0.)), + "Bad origin for default 2D line segment ctor."); + VTKM_MATH_ASSERT(test_equal(seg0.Endpoints[1], V2(1., 0.)), + "Bad direction for default 2D line segment ctor."); + + V2 p0(1., 1.); + V2 p1(3., 3.); + V2 p2(2., 2.); + // V2 p3(static_cast(1.2928932), static_cast(2.7071068)); + V2 dir(static_cast(-0.7071068), static_cast(0.7071068)); + vtkm::LineSegment seg1(p0, p1); + Line2 ray = seg1.PerpendicularBisector(); + VTKM_MATH_ASSERT(test_equal(ray.Origin, p2), "Perpendicular bisector origin failed in 2D."); + VTKM_MATH_ASSERT(test_equal(ray.Direction, dir), + "Perpendicular bisector direction failed in 2D."); + } + + { + using V3 = vtkm::Vec; + + vtkm::LineSegment seg0; + VTKM_MATH_ASSERT(test_equal(seg0.Endpoints[0], V3(0., 0., 0.)), + "Bad origin for default 3D line segment ctor."); + VTKM_MATH_ASSERT(test_equal(seg0.Endpoints[1], V3(1., 0., 0.)), + "Bad direction for default 3D line segment ctor."); + + V3 p0(1., 1., 0.); + V3 p1(3., 3., 0.); + V3 p2(2., 2., 0.); + V3 p3(static_cast(0.70710678), static_cast(0.70710678), 0.); + vtkm::LineSegment seg1(p0, p1); + vtkm::Plane bisector = seg1.PerpendicularBisector(); + VTKM_MATH_ASSERT(test_equal(bisector.Origin, p2), + "Perpendicular bisector origin failed in 3D."); + VTKM_MATH_ASSERT(test_equal(bisector.Normal, p3), + "Perpendicular bisector direction failed in 3D."); + } + + vtkm::Vec origin(0., 0., 0.); + for (vtkm::IdComponent index = 0; index < Coords::NUM_COORDS; ++index) + { + auto p0 = Coords{}.EndpointList(index); + auto p1 = Coords{}.EndpointList((index + 1) % Coords::NUM_COORDS); + + vtkm::LineSegment segment(p0, p1); + vtkm::Vec closest; + T param; + auto dp0 = segment.DistanceTo(p0); + auto dp1 = segment.DistanceTo(p1, param, closest); + VTKM_MATH_ASSERT(test_equal(dp0, 0.0), "Distance to endpoint 0 not zero."); + VTKM_MATH_ASSERT(test_equal(dp1, 0.0), "Distance to endpoint 1 not zero."); + VTKM_MATH_ASSERT(test_equal(param, 1.0), "Parameter value of endpoint 1 not 1.0."); + VTKM_MATH_ASSERT(test_equal(p1, closest), "Closest point not endpoint 1."); + + closest = segment.Evaluate(static_cast(0.0)); + VTKM_MATH_ASSERT(test_equal(p0, closest), "Evaluated point not endpoint 0."); + + auto dpo = segment.DistanceTo(origin, param, closest); + auto clo = Coords{}.ClosestToOriginList(index); + auto dst = Coords{}.DistanceToOriginList(index); + VTKM_MATH_ASSERT(test_equal(closest, clo), "Closest point to origin doesn't match."); + VTKM_MATH_ASSERT(test_equal(dpo, dst), "Distance to origin doesn't match."); + } + } +}; + +template +struct TryLineSegmentTests +{ + template + void operator()(const T&) const + { + vtkm::cont::DeviceAdapterAlgorithm::Schedule(LineSegmentTests(), 1); + } +}; + +//----------------------------------------------------------------------------- + +template +struct PlaneTests : public vtkm::exec::FunctorBase +{ + VTKM_EXEC + void operator()(vtkm::Id) const + { + vtkm::Vec origin(0., 0., 0.); + vtkm::Vec zvectr(0., 0., 5.); // intentionally not unit length to test normalization. + vtkm::Plane plane; + vtkm::LineSegment segment; + T dist; + bool didIntersect; + bool isLineInPlane; + vtkm::Vec nearest; + vtkm::Vec p0; + vtkm::Vec p1; + T param; + + // Test signed plane-point distance + plane = vtkm::Plane(origin, zvectr); + dist = plane.DistanceTo(vtkm::Vec(82., 0.5, 1.25)); + VTKM_MATH_ASSERT(test_equal(dist, 1.25), "Bad positive point-plane distance."); + dist = plane.DistanceTo(vtkm::Vec(82., 0.5, -1.25)); + VTKM_MATH_ASSERT(test_equal(dist, -1.25), "Bad negative point-plane distance."); + dist = plane.DistanceTo(vtkm::Vec(82., 0.5, 0.0)); + VTKM_MATH_ASSERT(test_equal(dist, 0.0), "Bad zero point-plane distance."); + + // Test line intersection + { + // Case 1. No intersection + segment = vtkm::LineSegment((p0 = vtkm::Vec(1., 1., 1.)), + (p1 = vtkm::Vec(2., 2., 2.))); + didIntersect = plane.Intersect(segment, param, nearest, isLineInPlane); + VTKM_MATH_ASSERT(test_equal(didIntersect, false), "Plane and line should not intersect (1)."); + VTKM_MATH_ASSERT(test_equal(isLineInPlane, false), + "Line improperly reported as in plane (1)."); + VTKM_MATH_ASSERT(test_equal(nearest, p0), "Unexpected nearest point (1)."); + VTKM_MATH_ASSERT(test_equal(param, 0.0), "Unexpected nearest parameter value (1)."); + + // Case 2. Degenerate intersection (entire segment lies in plane) + segment = vtkm::LineSegment((p0 = vtkm::Vec(1., 1., 0.)), + (p1 = vtkm::Vec(2., 2., 0.))); + didIntersect = plane.Intersect(segment, param, nearest, isLineInPlane); + VTKM_MATH_ASSERT(test_equal(didIntersect, true), "Plane and line should intersect (2)."); + VTKM_MATH_ASSERT(test_equal(isLineInPlane, true), + "Line improperly reported as out of plane (2)."); + + // Case 3. Endpoint intersection + segment = vtkm::LineSegment((p0 = vtkm::Vec(1., 1., 1.)), + (p1 = vtkm::Vec(2., 2., 0.))); + didIntersect = plane.Intersect(segment, param, nearest, isLineInPlane); + VTKM_MATH_ASSERT(test_equal(didIntersect, true), "Plane and line should intersect (3a)."); + VTKM_MATH_ASSERT(test_equal(isLineInPlane, false), + "Line improperly reported as in plane (3a)."); + VTKM_MATH_ASSERT(test_equal(param, 1.0), "Invalid parameter for intersection point (3a)."); + VTKM_MATH_ASSERT(test_equal(nearest, p1), "Invalid intersection point (3a)."); + + segment = vtkm::LineSegment((p0 = vtkm::Vec(1., 1., 0.)), + (p1 = vtkm::Vec(2., 2., 1.))); + didIntersect = plane.Intersect(segment, param, nearest, isLineInPlane); + VTKM_MATH_ASSERT(test_equal(didIntersect, true), "Plane and line should intersect (3b)."); + VTKM_MATH_ASSERT(test_equal(isLineInPlane, false), + "Line improperly reported as in plane (3b)."); + VTKM_MATH_ASSERT(test_equal(param, 0.0), "Invalid parameter for intersection point (3b)."); + VTKM_MATH_ASSERT(test_equal(nearest, p0), "Invalid intersection point (3b)."); + + // Case 4. General-position intersection + segment = vtkm::LineSegment((p0 = vtkm::Vec(-1., -1., -1.)), + (p1 = vtkm::Vec(2., 2., 1.))); + didIntersect = plane.Intersect(segment, param, nearest, isLineInPlane); + VTKM_MATH_ASSERT(test_equal(didIntersect, true), "Plane and line should intersect (4)."); + VTKM_MATH_ASSERT(test_equal(isLineInPlane, false), + "Line improperly reported as in plane (4)."); + VTKM_MATH_ASSERT(test_equal(param, 0.5), "Invalid parameter for intersection point (4)."); + VTKM_MATH_ASSERT(test_equal(nearest, vtkm::Vec(0.5, 0.5, 0)), + "Invalid intersection point (4)."); + } + + // Test plane-plane intersection + { + using V3 = vtkm::Vec; + using PlaneType = vtkm::Plane; + // Case 1. Coincident planes + p0 = V3(1., 2., 3.); + p1 = V3(5., 7., -6.); + V3 nn = vtkm::Normal(V3(1., 1., 1)); + PlaneType pa(p0, nn); + PlaneType pb(p1, nn); + vtkm::Line3 ii; + bool coincident; + didIntersect = pa.Intersect(pb, ii, coincident); + VTKM_MATH_ASSERT(test_equal(didIntersect, false), + "Coincident planes should have degenerate intersection."); + VTKM_MATH_ASSERT(test_equal(coincident, true), + "Coincident planes should be marked coincident."); + + // Case 2. Offset planes + p1 = V3(5., 6., 7.); + pb = PlaneType(p1, nn); + didIntersect = pa.Intersect(pb, ii, coincident); + VTKM_MATH_ASSERT(test_equal(didIntersect, false), + "Offset planes should have degenerate intersection."); + VTKM_MATH_ASSERT(test_equal(coincident, false), + "Offset planes should not be marked coincident."); + + // Case 3. General position + p1 = V3(1., 2., 0.); + V3 n2(0., 0., 1.); + pb = PlaneType(p1, n2); + didIntersect = pa.Intersect(pb, ii, coincident); + VTKM_MATH_ASSERT(test_equal(didIntersect, true), + "Proper planes should have non-degenerate intersection."); + VTKM_MATH_ASSERT(test_equal(coincident, false), + "Proper planes should not be marked coincident."); + VTKM_MATH_ASSERT(test_equal(ii.Origin, V3(2.5, 3.5, 0)), + "Unexpected intersection-line base point."); + VTKM_MATH_ASSERT(test_equal(ii.Direction, vtkm::Normal(V3(1, -1, 0))), + "Unexpected intersection-line direction."); + } + } +}; + +template +struct TryPlaneTests +{ + template + void operator()(const T&) const + { + vtkm::cont::DeviceAdapterAlgorithm::Schedule(PlaneTests(), 1); + } +}; + +//----------------------------------------------------------------------------- + +template +struct SphereTests : public vtkm::exec::FunctorBase +{ + VTKM_EXEC + void operator()(vtkm::Id) const + { + { + using V2 = vtkm::Vec; + V2 origin(0., 0.); + vtkm::Sphere defaultSphere; + VTKM_MATH_ASSERT(test_equal(defaultSphere.Center, origin), "Default circle not at origin."); + VTKM_MATH_ASSERT(test_equal(defaultSphere.Radius, 1.0), "Default circle not unit radius."); + + vtkm::Sphere sphere(origin, -2.); + VTKM_MATH_ASSERT(test_equal(sphere.Radius, -1.0), "Negative radius should be reset to -1."); + VTKM_MATH_ASSERT(test_equal(sphere.IsValid(), false), + "Negative radius should leave sphere invalid."); + + sphere = vtkm::Circle(origin, 1.0); + VTKM_MATH_ASSERT(test_equal(sphere.IsValid(), true), "Circle assignment failed."); + VTKM_MATH_ASSERT(test_equal(sphere.Contains(origin), true), + "Circle does not contain its center."); + VTKM_MATH_ASSERT(test_equal(sphere.Classify(V2(1., 0.)), 0), "Circle point not on boundary."); + VTKM_MATH_ASSERT(test_equal(sphere.Classify(V2(0.75, 0.75)), +1), + "Circle contains a point that should be outside."); + + V2 p0(static_cast(-0.7071), static_cast(-0.7071)); + V2 p1(static_cast(+0.7071), static_cast(-0.7071)); + V2 p2(static_cast(0.0), static_cast(1.0)); + sphere = make_CircleFrom3Points(p0, p1, p2); + VTKM_MATH_ASSERT(test_equal(sphere.IsValid(), true), "Could not create 3-point circle."); + + V2 p3(1, 1); + V2 p4(3, 4); + V2 p5(5, 12); + sphere = make_CircleFrom3Points(p3, p4, p5); + VTKM_MATH_ASSERT(test_equal(sphere.IsValid(), true), "Could not create 3-point circle."); + T tol = static_cast(1e-3); // Use a loose tolerance + VTKM_MATH_ASSERT(test_equal(sphere.Center, vtkm::Vec(-12.4f, 12.1f)), + "Invalid circle center."); + VTKM_MATH_ASSERT(test_equal(sphere.Radius, static_cast(17.400291f)), + "Invalid circle radius."); + VTKM_MATH_ASSERT(test_equal(sphere.Classify(p3, tol), 0), + "Generator p3 not on circle boundary."); + VTKM_MATH_ASSERT(test_equal(sphere.Classify(p4, tol), 0), + "Generator p4 not on circle boundary."); + VTKM_MATH_ASSERT(test_equal(sphere.Classify(p5, tol), 0), + "Generator p5 not on circle boundary."); + + V2 p6(1, 1); + V2 p7(4, 4); + V2 p8(5, 5); + sphere = make_CircleFrom3Points(p6, p7, p8); + VTKM_MATH_ASSERT(test_equal(sphere.IsValid(), false), + "3-point circle construction should fail with points on a line."); + } + { + using V3 = vtkm::Vec; + + V3 p0(0., 1., 0.); + V3 p1(1., 0., 0.); + V3 p2(-1., 0., 0.); + V3 p3(0., 0., 1.); + V3 p4 = vtkm::Normal(V3(1., 1., 1.)); + + V3 origin(0., 0., 0.); + vtkm::Sphere defaultSphere; + VTKM_MATH_ASSERT(test_equal(defaultSphere.Center, origin), "Default sphere not at origin."); + VTKM_MATH_ASSERT(test_equal(defaultSphere.Radius, 1.0), "Default sphere not unit radius."); + + vtkm::Sphere sphere = make_SphereFrom4Points(p0, p1, p2, p3, static_cast(1.0e-6)); + VTKM_MATH_ASSERT(test_equal(sphere.IsValid(), true), "Easy sphere 1 not valid."); + VTKM_MATH_ASSERT(test_equal(sphere.Center, origin), "Easy sphere 1 not at origin."); + VTKM_MATH_ASSERT(test_equal(sphere.Radius, 1.0), "Easy sphere 1 not unit radius."); + + sphere = make_SphereFrom4Points(p0, p1, p2, p4, static_cast(1.0e-6)); + VTKM_MATH_ASSERT(test_equal(sphere.IsValid(), true), "Easy sphere 2 not valid."); + VTKM_MATH_ASSERT(test_equal(sphere.Center, origin), "Easy sphere 2 not at origin."); + VTKM_MATH_ASSERT(test_equal(sphere.Radius, 1.0), "Easy sphere 2 not unit radius."); + + V3 fancyCenter(1, 2, 3); + T fancyRadius(2.5); + + V3 fp0 = fancyCenter + fancyRadius * p0; + V3 fp1 = fancyCenter + fancyRadius * p1; + V3 fp2 = fancyCenter + fancyRadius * p2; + V3 fp4 = fancyCenter + fancyRadius * p4; + + sphere = make_SphereFrom4Points(fp0, fp1, fp2, fp4, static_cast(1.0e-6)); + VTKM_MATH_ASSERT(test_equal(sphere.IsValid(), true), "Medium sphere 1 not valid."); + VTKM_MATH_ASSERT(test_equal(sphere.Center, fancyCenter), "Medium sphere 1 not at (1,2,3)."); + VTKM_MATH_ASSERT(test_equal(sphere.Radius, fancyRadius), "Medium sphere 1 not radius 2.5."); + } + } +}; + +template +struct TrySphereTests +{ + template + void operator()(const T&) const + { + vtkm::cont::DeviceAdapterAlgorithm::Schedule(SphereTests(), 1); + } +}; + +//----------------------------------------------------------------------------- +template +void RunGeometryTests() +{ + std::cout << "Tests for rays." << std::endl; + vtkm::testing::Testing::TryTypes(TryRayTests(), vtkm::TypeListTagFieldScalar()); + std::cout << "Tests for line segments." << std::endl; + vtkm::testing::Testing::TryTypes(TryLineSegmentTests(), vtkm::TypeListTagFieldScalar()); + std::cout << "Tests for planes." << std::endl; + vtkm::testing::Testing::TryTypes(TryPlaneTests(), vtkm::TypeListTagFieldScalar()); + std::cout << "Tests for spheres." << std::endl; + vtkm::testing::Testing::TryTypes(TrySphereTests(), vtkm::TypeListTagFieldScalar()); +} + +} // namespace UnitTestGeometryNamespace + +#endif //vtk_m_testing_TestingGeometry_h diff --git a/vtkm/testing/UnitTestVectorAnalysis.cxx b/vtkm/testing/UnitTestVectorAnalysis.cxx index 4cceb93c3..c356c1773 100644 --- a/vtkm/testing/UnitTestVectorAnalysis.cxx +++ b/vtkm/testing/UnitTestVectorAnalysis.cxx @@ -158,6 +158,21 @@ void TestCross(const vtkm::Vec& x, const vtkm::Vec& y) "Triangle normal is not really normal."); } +template +void TestOrthonormalize(const VectorBasisType& inputs, int expectedRank) +{ + VectorBasisType outputs; + int actualRank = vtkm::Orthonormalize(inputs, outputs); + std::cout << "Testing orthonormalize\n" + << " Rank " << actualRank << " expected " << expectedRank << "\n" + << " Basis vectors:\n"; + for (int i = 0; i < actualRank; ++i) + { + std::cout << " " << i << " " << outputs[i] << "\n"; + } + VTKM_TEST_ASSERT(test_equal(actualRank, expectedRank), "Orthonormalized rank is unexpected."); +} + struct TestLinearFunctor { template @@ -208,10 +223,43 @@ struct TestCrossFunctor } }; +struct TestVectorFunctor +{ + template + void operator()(const VectorType&) const + { + constexpr int NUM_COMPONENTS = VectorType::NUM_COMPONENTS; + using Traits = vtkm::VecTraits; + using ComponentType = typename Traits::ComponentType; + vtkm::Vec basis; + VectorType normalizedVector = VectorType(vtkm::RSqrt(ComponentType(NUM_COMPONENTS))); + VectorType zeroVector = VectorType(ComponentType(0)); + // Test with a degenerate set of inputs: + basis[0] = zeroVector; + basis[1] = normalizedVector; + for (int ii = 2; ii < NUM_COMPONENTS; ++ii) + { + basis[ii] = zeroVector; + } + TestOrthonormalize(basis, 1); + // Now test with a valid set of inputs: + for (int ii = 0; ii < NUM_COMPONENTS; ++ii) + { + for (int jj = 0; jj < NUM_COMPONENTS; ++jj) + { + basis[ii][jj] = + ComponentType(jj == ii ? 1.0 : 0.0) + ComponentType(0.05) * ComponentType(jj); + } + } + TestOrthonormalize(basis, NUM_COMPONENTS); + } +}; + void TestVectorAnalysis() { vtkm::testing::Testing::TryTypes(TestLinearFunctor(), vtkm::TypeListTagField()); vtkm::testing::Testing::TryTypes(TestCrossFunctor(), vtkm::TypeListTagFieldVec3()); + vtkm::testing::Testing::TryTypes(TestVectorFunctor(), vtkm::TypeListTagFloatVec()); } } // anonymous namespace