From e909edd61ddf5b96027619cc7a0d2ac25035b7b1 Mon Sep 17 00:00:00 2001 From: VTK-c Upstream Date: Tue, 24 Sep 2019 21:20:15 -0400 Subject: [PATCH] vtkc 2019-09-24 (cdc72582) Code extracted from: https://gitlab.kitware.com/sujin.philip/vtk-c.git at commit cdc72582dae24cc66b3525961e32edbc2a395341 (master). --- LICENSE.md | 43 +++ README.md | 65 ++++ vtkc/ErrorCode.h | 69 ++++ vtkc/FieldAccessor.h | 367 +++++++++++++++++++++ vtkc/Hexahedron.h | 254 +++++++++++++++ vtkc/Line.h | 176 ++++++++++ vtkc/Pixel.h | 178 +++++++++++ vtkc/Polygon.h | 541 +++++++++++++++++++++++++++++++ vtkc/Pyramid.h | 293 +++++++++++++++++ vtkc/Quad.h | 195 +++++++++++ vtkc/Shapes.h | 142 ++++++++ vtkc/Tetra.h | 208 ++++++++++++ vtkc/Triangle.h | 271 ++++++++++++++++ vtkc/Vertex.h | 131 ++++++++ vtkc/Voxel.h | 136 ++++++++ vtkc/Wedge.h | 228 +++++++++++++ vtkc/internal/Common.h | 310 ++++++++++++++++++ vtkc/internal/Config.h | 38 +++ vtkc/internal/Math.h | 712 +++++++++++++++++++++++++++++++++++++++++ vtkc/vtkc.h | 251 +++++++++++++++ 20 files changed, 4608 insertions(+) create mode 100644 LICENSE.md create mode 100644 README.md create mode 100644 vtkc/ErrorCode.h create mode 100644 vtkc/FieldAccessor.h create mode 100644 vtkc/Hexahedron.h create mode 100644 vtkc/Line.h create mode 100644 vtkc/Pixel.h create mode 100644 vtkc/Polygon.h create mode 100644 vtkc/Pyramid.h create mode 100644 vtkc/Quad.h create mode 100644 vtkc/Shapes.h create mode 100644 vtkc/Tetra.h create mode 100644 vtkc/Triangle.h create mode 100644 vtkc/Vertex.h create mode 100644 vtkc/Voxel.h create mode 100644 vtkc/Wedge.h create mode 100644 vtkc/internal/Common.h create mode 100644 vtkc/internal/Config.h create mode 100644 vtkc/internal/Math.h create mode 100644 vtkc/vtkc.h diff --git a/LICENSE.md b/LICENSE.md new file mode 100644 index 000000000..c9bb454d8 --- /dev/null +++ b/LICENSE.md @@ -0,0 +1,43 @@ +# VTK-c License Version 1.0 # + +Copyright (c) 2019 +Kitware Inc., +National Technology & Engineering Solutions of Sandia, LLC (NTESS), +UT-Battelle, LLC., +Los Alamos National Security, LLC., +All rights reserved. + +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. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the + distribution. + + * Neither the name of Kitware nor the names of any contributors may + be used to endorse or promote products derived from this software + without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/README.md b/README.md new file mode 100644 index 000000000..e2231d7db --- /dev/null +++ b/README.md @@ -0,0 +1,65 @@ +# VTK-c # + +VTK-c is a lightweight collection of cell types and cell functionality +that is designed be used scientific visualization libraries of +VTK-m, and VTK. + +You can find out more about the design of VTK-c in [DESIGN.md]. + +## Contributing ## + +There are many ways to contribute to [VTK-c]: + + + Submit new or add to discussions of a feature requests or bugs on the + [VTK-c Issue Tracker]. + + + Submit a Pull Request + + See [CONTRIBUTING.md] for detailed instructions on how to create a + Pull Request. + + See the [VTK-c Coding Conventions] that must be followed for + contributed code. + +## Compiler Requirements ## + + + C++11 Compiler. VTK-c has been confirmed to work with the following + + GCC 4.8+ + + Clang 3.3+ + + XCode 5.0+ + + MSVC 2015+ + +## Example## + +Below is a simple example of using VTK-c to get derivatives and +parametric coordinates for different cell types: + +```cpp +#incude + +std::array pcoords; +auto status = vtkc::parametricCenter(vtkc::Hexahedron{}, pcoords); + +std::vector> points = { ... }; +std::vector> field = { ... }; +std::array derivs[3]; +status = vtkc::derivative(vtkc::Hexahedron{}, + vtkc::makeFieldAccessorNestedSOAConst(points, 3), + vtkc::makeFieldAccessorNestedSOAConst(field, 4), + pcoords, + derivs[0], + derivs[1], + derivs[2]); +``` + +## License ## + +VTK-c is distributed under the OSI-approved BSD 3-clause License. +See [LICENSE.md] for details. + + +[VTK-c]: https://gitlab.kitware.com/sujin.philip/vtk-c/ +[VTK-c Issue Tracker]: https://gitlab.kitware.com/sujin.philip/vtk-c/issues + +[CONTRIBUTING.md]: CONTRIBUTING.md +[DESIGN.md]: docs/Design.md +[LICENSE.md]: LICENSE.md +[VTK-c Coding Conventions]: docs/CodingConventions.md diff --git a/vtkc/ErrorCode.h b/vtkc/ErrorCode.h new file mode 100644 index 000000000..eedfcbb33 --- /dev/null +++ b/vtkc/ErrorCode.h @@ -0,0 +1,69 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.md 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. +//============================================================================ +#ifndef vtk_c_ErrorCode_h +#define vtk_c_ErrorCode_h + +#include + +#include + +namespace vtkc +{ + +enum class ErrorCode : std::int32_t +{ + SUCCESS = 0, + INVALID_SHAPE_ID, + INVALID_NUMBER_OF_POINTS, + WRONG_SHAPE_ID_FOR_TAG_TYPE, + INVALID_POINT_ID, + SOLUTION_DID_NOT_CONVERGE, + MATRIX_LUP_FACTORIZATION_FAILED, + DEGENERATE_CELL_DETECTED +}; + +VTKC_EXEC +inline const char* errorString(ErrorCode code) noexcept +{ + switch (code) + { + case ErrorCode::SUCCESS: + return "Success"; + case ErrorCode::INVALID_SHAPE_ID: + return "Invalid shape id"; + case ErrorCode::INVALID_NUMBER_OF_POINTS: + return "Invalid number of points"; + case ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE: + return "Wrong shape id for tag type"; + case ErrorCode::INVALID_POINT_ID: + return "Invalid point id"; + case ErrorCode::SOLUTION_DID_NOT_CONVERGE: + return "Solution did not converge"; + case ErrorCode::MATRIX_LUP_FACTORIZATION_FAILED: + return "LUP factorization failed"; + case ErrorCode::DEGENERATE_CELL_DETECTED: + return "Degenerate cell detected"; + } + + return "Invalid error"; +} + +} // vtkc + +#define VTKC_RETURN_ON_ERROR(call) \ + { \ + auto status = call; \ + if (status != vtkc::ErrorCode::SUCCESS) \ + { \ + return status; \ + } \ + } + +#endif // vtk_c_ErrorCode_h diff --git a/vtkc/FieldAccessor.h b/vtkc/FieldAccessor.h new file mode 100644 index 000000000..d6252be80 --- /dev/null +++ b/vtkc/FieldAccessor.h @@ -0,0 +1,367 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.md 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. +//============================================================================ +#ifndef vtk_c_FieldAccessor_h +#define vtk_c_FieldAccessor_h + +#include + +#include + +namespace vtkc +{ +namespace internal +{ +//---------------------------------------------------------------------------- +/// Does T support the indexing operator? +/// Examples: +/// std::array +/// std::vector +/// T* +/// custom classes with operator[](integer_type) +/// + +struct NotVecTag {}; + +template inline decltype(std::declval()[0]) VecCheck(T); +inline NotVecTag VecCheck(...); + +template +using IsVecType = + std::integral_constant())), + internal::NotVecTag>::value>; + +//---------------------------------------------------------------------------- +/// Does T support operator()(integer_type)? + +struct NotFuncTag {}; + +template inline decltype(std::declval()(0)) FuncCheck(T); +inline NotFuncTag FuncCheck(...); + +template +using IsFuncType = + std::integral_constant())), + internal::NotFuncTag>::value>; + +//---------------------------------------------------------------------------- +template typename T::size_type IndexTypeImpl(T); +IdComponent IndexTypeImpl(...); + +template +using IndexType = decltype(IndexTypeImpl(std::declval())); + +//---------------------------------------------------------------------------- +struct VecTypeIndexer {}; +struct VecTypeFunctor {}; +struct VecTypeScalar {}; + +template +using VecTypeTag = + typename std::conditional::value, + VecTypeIndexer, + typename std::conditional::value, + VecTypeFunctor, + VecTypeScalar>::type>::type; + +//---------------------------------------------------------------------------- +template +VTKC_EXEC auto componentImpl(VecTypeIndexer, T&& vec, int idx) -> decltype(vec[0]) +{ + return vec[static_cast>(idx)]; +} + +template +VTKC_EXEC auto componentImpl(VecTypeFunctor, T&& vec, int idx) -> decltype(vec(0)) +{ + return vec(idx); +} + +template +VTKC_EXEC T& componentImpl(VecTypeScalar, T&& vec, int) +{ + return vec; +} + +} // namespace internal + +//---------------------------------------------------------------------------- +template +VTKC_EXEC auto component(T&& vec, int idx) + -> decltype(internal::componentImpl(internal::VecTypeTag{}, std::forward(vec), 0)) +{ + return internal::componentImpl(internal::VecTypeTag{}, std::forward(vec), idx); +} + +//---------------------------------------------------------------------------- +template +using ComponentType = typename std::decay(), 0))>::type; + +///============================================================================ +/// Since there are different ways fields maybe represented in the clients of +/// this libarary, VTK-c relies on helper classes that implement the +/// \c FieldAccessor "concept" to access the elements of a field. +/// +/// These classes should wrap the field and provide the follolwing interface: +/// +/// template +/// class FieldAccessor +/// { +/// public: +/// /// An alias for the component type of the field +/// using ValueType = ...; +/// +/// /// Return the number of components +/// int getNumberOfComponents() const; +/// +/// /// Set the value to `tuple` and `comp` to `value`. +/// void setValue(int tuple, int comp, const ValueType& value) const; +/// +/// /// Get the value at `tuple` and `comp`. +/// ValueType getValue(int tuple, int comp) const; +/// +/// /// Set the tuple at index `tuple`. It is recomended to make this a +/// /// template function and use `vtkc::component` to access +/// /// the components of `value`. +/// template +/// void setTuple(int tuple, const VecType& value) const; +/// +/// /// Get the tuple at index `tuple`. It is recomended to make this a +/// /// template function and use `vtkc::component` to access +/// /// the components of `value`. +/// template +/// void getTuple(int tuple, VecType& value) const; +/// } +/// +/// The set functions are optional and such a class would act as a const variant +/// of \c FieldAccessor. +/// + +///---------------------------------------------------------------------------- +template +class FieldAccessorNestedSOA +{ +public: + using ValueType = typename std::decay()[0], 0))>::type; + + VTKC_EXEC FieldAccessorNestedSOA(FieldType& field, int numberOfComponents = 1) + : Field(&field), NumberOfComponents(numberOfComponents) + { + } + + VTKC_EXEC int getNumberOfComponents() const + { + return this->NumberOfComponents; + } + + VTKC_EXEC void setValue(int tuple, int comp, const ValueType& value) const + { + component((*this->Field)[static_cast(tuple)], comp) = value; + } + + VTKC_EXEC ValueType getValue(int tuple, int comp) const + { + return component((*this->Field)[static_cast(tuple)], comp); + } + + template + VTKC_EXEC void setTuple(int tuple, const VecType& value) const + { + for (int i = 0; i < this->NumberOfComponents; ++i) + { + component((*this->Field)[static_cast(tuple)], i) = + static_cast(component(value, i)); + } + } + + template + VTKC_EXEC void getTuple(int tuple, VecType& value) const + { + for (int i = 0; i < this->NumberOfComponents; ++i) + { + component(value, i) = static_cast>( + component((*this->Field)[static_cast(tuple)], i)); + } + } + +private: + using IdxType = internal::IndexType; + + FieldType* Field; + int NumberOfComponents; +}; + +template +class FieldAccessorNestedSOAConst +{ +public: + using ValueType = typename std::decay()[0], 0))>::type; + + VTKC_EXEC FieldAccessorNestedSOAConst(const FieldType& field, int numberOfComponents = 1) + : Field(&field), NumberOfComponents(numberOfComponents) + { + } + + VTKC_EXEC int getNumberOfComponents() const + { + return this->NumberOfComponents; + } + + VTKC_EXEC ValueType getValue(int tuple, int comp) const + { + return component((*this->Field)[static_cast(tuple)], comp); + } + + template + VTKC_EXEC void getTuple(int tuple, VecType& value) const + { + for (int i = 0; i < this->NumberOfComponents; ++i) + { + component(value, i) = static_cast>( + component((*this->Field)[static_cast(tuple)], i)); + } + } + +private: + using IdxType = internal::IndexType; + + const FieldType* Field; + int NumberOfComponents; +}; + +template +VTKC_EXEC +FieldAccessorNestedSOA makeFieldAccessorNestedSOA(FieldType& field, + int numberOfComponents = 1) +{ + return FieldAccessorNestedSOA(field, numberOfComponents); +} + +template +VTKC_EXEC +FieldAccessorNestedSOAConst makeFieldAccessorNestedSOAConst(const FieldType& field, + int numberOfComponents = 1) +{ + return FieldAccessorNestedSOAConst(field, numberOfComponents); +} + +///---------------------------------------------------------------------------- +template +class FieldAccessorFlatSOA +{ +public: + using ValueType = typename std::decay()[0])>::type; + + VTKC_EXEC FieldAccessorFlatSOA(FieldType& field, int numberOfComponents = 1) + : Field(&field), NumberOfComponents(numberOfComponents) + { + } + + VTKC_EXEC int getNumberOfComponents() const + { + return this->NumberOfComponents; + } + + VTKC_EXEC void setValue(int tuple, int comp, const ValueType& value) const + { + auto FlatIdx = static_cast(tuple * this->NumberOfComponents + comp); + (*this->Field)[FlatIdx] = value; + } + + VTKC_EXEC ValueType getValue(int tuple, int comp) const + { + auto FlatIdx = static_cast(tuple * this->NumberOfComponents + comp); + return (*this->Field)[FlatIdx]; + } + + template + VTKC_EXEC void setTuple(int tuple, const VecType& value) const + { + auto start = static_cast(tuple * this->NumberOfComponents); + for (int i = 0; i < this->NumberOfComponents; ++i) + { + (*this->Field)[start++] = static_cast(component(value, i)); + } + } + + template + VTKC_EXEC void getTuple(int tuple, VecType& value) const + { + auto start = static_cast(tuple * this->NumberOfComponents); + for (int i = 0; i < this->NumberOfComponents; ++i) + { + component(value, i) = static_cast>((*this->Field)[start++]); + } + } + +private: + using IdxType = internal::IndexType; + + FieldType* Field; + int NumberOfComponents; +}; + +template +class FieldAccessorFlatSOAConst +{ +public: + using ValueType = typename std::decay()[0])>::type; + + VTKC_EXEC FieldAccessorFlatSOAConst(const FieldType& field, int numberOfComponents = 1) + : Field(&field), NumberOfComponents(numberOfComponents) + { + } + + VTKC_EXEC int getNumberOfComponents() const + { + return this->NumberOfComponents; + } + + VTKC_EXEC ValueType getValue(int tuple, int comp) const + { + auto FlatIdx = static_cast(tuple * this->NumberOfComponents + comp); + return (*this->Field)[FlatIdx]; + } + + template + VTKC_EXEC void getTuple(int tuple, VecType& value) const + { + auto start = static_cast(tuple * this->NumberOfComponents); + for (int i = 0; i < this->NumberOfComponents; ++i) + { + component(value, i) = static_cast>((*this->Field)[start++]); + } + } + +private: + using IdxType = internal::IndexType; + + const FieldType* Field; + int NumberOfComponents; +}; + +template +VTKC_EXEC +FieldAccessorFlatSOA makeFieldAccessorFlatSOA(FieldType& field, + int numberOfComponents) +{ + return FieldAccessorFlatSOA(field, numberOfComponents); +} + +template +VTKC_EXEC +FieldAccessorFlatSOAConst makeFieldAccessorFlatSOAConst(const FieldType& field, + int numberOfComponents) +{ + return FieldAccessorFlatSOAConst(field, numberOfComponents); +} + +} // namespace vtkc + +#endif // vtk_c_FieldAccessor_h diff --git a/vtkc/Hexahedron.h b/vtkc/Hexahedron.h new file mode 100644 index 000000000..6717d5fd6 --- /dev/null +++ b/vtkc/Hexahedron.h @@ -0,0 +1,254 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.md 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. +//============================================================================ +#ifndef vtk_c_Hexahedron_h +#define vtk_c_Hexahedron_h + +#include +#include + +#include + +namespace vtkc +{ + +class Hexahedron : public Cell +{ +public: + constexpr VTKC_EXEC Hexahedron() : Cell(ShapeId::HEXAHEDRON, 8) {} + constexpr VTKC_EXEC explicit Hexahedron(const Cell& cell) : Cell(cell) {} +}; + +VTKC_EXEC inline vtkc::ErrorCode validate(Hexahedron tag) noexcept +{ + if (tag.shape() != ShapeId::HEXAHEDRON && tag.shape() != ShapeId::VOXEL) + { + return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE; + } + if (tag.numberOfPoints() != 8) + { + return ErrorCode::INVALID_NUMBER_OF_POINTS; + } + + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricCenter(Hexahedron, CoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + component(pcoords, 0) = 0.5f; + component(pcoords, 1) = 0.5f; + component(pcoords, 2) = 0.5f; + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricPoint( + Hexahedron, IdComponent pointId, CoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + switch (pointId) + { + case 0: + component(pcoords, 0) = 0.0f; + component(pcoords, 1) = 0.0f; + component(pcoords, 2) = 0.0f; + break; + case 1: + component(pcoords, 0) = 1.0f; + component(pcoords, 1) = 0.0f; + component(pcoords, 2) = 0.0f; + break; + case 2: + component(pcoords, 0) = 1.0f; + component(pcoords, 1) = 1.0f; + component(pcoords, 2) = 0.0f; + break; + case 3: + component(pcoords, 0) = 0.0f; + component(pcoords, 1) = 1.0f; + component(pcoords, 2) = 0.0f; + break; + case 4: + component(pcoords, 0) = 0.0f; + component(pcoords, 1) = 0.0f; + component(pcoords, 2) = 1.0f; + break; + case 5: + component(pcoords, 0) = 1.0f; + component(pcoords, 1) = 0.0f; + component(pcoords, 2) = 1.0f; + break; + case 6: + component(pcoords, 0) = 1.0f; + component(pcoords, 1) = 1.0f; + component(pcoords, 2) = 1.0f; + break; + case 7: + component(pcoords, 0) = 0.0f; + component(pcoords, 1) = 1.0f; + component(pcoords, 2) = 1.0f; + break; + default: + return ErrorCode::INVALID_POINT_ID; + } + + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline ComponentType parametricDistance(Hexahedron, const CoordType& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + return internal::findParametricDistance(pcoords, 3); +} + +template +VTKC_EXEC inline bool cellInside(Hexahedron, const CoordType& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + using T = ComponentType; + return component(pcoords, 0) >= T{0} && component(pcoords, 0) <= T{1} && + component(pcoords, 1) >= T{0} && component(pcoords, 1) <= T{1} && + component(pcoords, 2) >= T{0} && component(pcoords, 2) <= T{1}; +} + +template +VTKC_EXEC inline vtkc::ErrorCode interpolate( + Hexahedron, + const Values& values, + const CoordType& pcoords, + Result&& result) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + using T = internal::ClosestFloatType; + + for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c) + { + auto vbf = internal::lerp(static_cast(values.getValue(0, c)), + static_cast(values.getValue(1, c)), + static_cast(component(pcoords, 0))); + auto vbb = internal::lerp(static_cast(values.getValue(3, c)), + static_cast(values.getValue(2, c)), + static_cast(component(pcoords, 0))); + auto vtf = internal::lerp(static_cast(values.getValue(4, c)), + static_cast(values.getValue(5, c)), + static_cast(component(pcoords, 0))); + auto vtb = internal::lerp(static_cast(values.getValue(7, c)), + static_cast(values.getValue(6, c)), + static_cast(component(pcoords, 0))); + auto vb = internal::lerp(vbf, vbb, static_cast(component(pcoords, 1))); + auto vt = internal::lerp(vtf, vtb, static_cast(component(pcoords, 1))); + auto v = internal::lerp(vb, vt, static_cast(component(pcoords, 2))); + component(result, c) = static_cast>(v); + } + + return ErrorCode::SUCCESS; +} + +namespace internal +{ + +template +VTKC_EXEC inline void parametricDerivative(Hexahedron, + const Values& values, + IdComponent comp, + const CoordType& pcoords, + Result&& result) noexcept +{ + using T = internal::ClosestFloatType; + T p0 = static_cast(component(pcoords, 0)); + T p1 = static_cast(component(pcoords, 1)); + T p2 = static_cast(component(pcoords, 2)); + T rm = T{1} - p0; + T sm = T{1} - p1; + T tm = T{1} - p2; + + T dr = (static_cast(values.getValue(0, comp)) * -sm * tm) + + (static_cast(values.getValue(1, comp)) * sm * tm) + + (static_cast(values.getValue(2, comp)) * p1 * tm) + + (static_cast(values.getValue(3, comp)) * -p1 * tm) + + (static_cast(values.getValue(4, comp)) * -sm * p2) + + (static_cast(values.getValue(5, comp)) * sm * p2) + + (static_cast(values.getValue(6, comp)) * p1 * p2) + + (static_cast(values.getValue(7, comp)) * -p1 * p2); + + T ds = (static_cast(values.getValue(0, comp)) * -rm * tm) + + (static_cast(values.getValue(1, comp)) * -p0 * tm) + + (static_cast(values.getValue(2, comp)) * p0 * tm) + + (static_cast(values.getValue(3, comp)) * rm * tm) + + (static_cast(values.getValue(4, comp)) * -rm * p2) + + (static_cast(values.getValue(5, comp)) * -p0 * p2) + + (static_cast(values.getValue(6, comp)) * p0 * p2) + + (static_cast(values.getValue(7, comp)) * rm * p2); + + T dt = (static_cast(values.getValue(0, comp)) * -rm * sm) + + (static_cast(values.getValue(1, comp)) * -p0 * sm) + + (static_cast(values.getValue(2, comp)) * -p0 * p1) + + (static_cast(values.getValue(3, comp)) * -rm * p1) + + (static_cast(values.getValue(4, comp)) * rm * sm) + + (static_cast(values.getValue(5, comp)) * p0 * sm) + + (static_cast(values.getValue(6, comp)) * p0 * p1) + + (static_cast(values.getValue(7, comp)) * rm * p1); + + component(result, 0) = static_cast>(dr); + component(result, 1) = static_cast>(ds); + component(result, 2) = static_cast>(dt); +} + +} // internal + +template +VTKC_EXEC inline vtkc::ErrorCode derivative( + Hexahedron, + const Points& points, + const Values& values, + const CoordType& pcoords, + Result&& dx, + Result&& dy, + Result&& dz) noexcept +{ + return internal::derivative3D(Hexahedron{}, + points, + values, + pcoords, + std::forward(dx), + std::forward(dy), + std::forward(dz)); +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricToWorld( + Hexahedron, + const Points& points, + const PCoordType& pcoords, + WCoordType&& wcoords) noexcept +{ + return interpolate(Hexahedron{}, points, pcoords, std::forward(wcoords)); +} + +template +VTKC_EXEC inline vtkc::ErrorCode worldToParametric( + Hexahedron, + const Points& points, + const WCoordType& wcoords, + PCoordType&& pcoords) noexcept +{ + return internal::worldToParametric3D( + Hexahedron{}, points, wcoords, std::forward(pcoords)); +} + +} // vtkc + +#endif // vtk_c_Hexahedron_h diff --git a/vtkc/Line.h b/vtkc/Line.h new file mode 100644 index 000000000..aa5134139 --- /dev/null +++ b/vtkc/Line.h @@ -0,0 +1,176 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.md 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. +//============================================================================ +#ifndef vtk_c_Line_h +#define vtk_c_Line_h + +#include +#include + +#include + +namespace vtkc +{ + +class Line : public Cell +{ +public: + constexpr VTKC_EXEC Line() : Cell(ShapeId::LINE, 2) {} + constexpr VTKC_EXEC explicit Line(const Cell& cell) : Cell(cell) {} +}; + +VTKC_EXEC inline vtkc::ErrorCode validate(Line tag) noexcept +{ + if (tag.shape() != ShapeId::LINE) + { + return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE; + } + if (tag.numberOfPoints() != 2) + { + return ErrorCode::INVALID_NUMBER_OF_POINTS; + } + + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricCenter(Line, CoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + component(pcoords, 0) = 0.5f; + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricPoint( + Line, IdComponent pointId, CoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + switch (pointId) + { + case 0: + component(pcoords, 0) = 0.0f; + return ErrorCode::SUCCESS; + case 1: + component(pcoords, 0) = 1.0f; + return ErrorCode::SUCCESS; + default: + return ErrorCode::INVALID_POINT_ID; + } +} + +template +VTKC_EXEC inline ComponentType parametricDistance(Line, const CoordType& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + return internal::findParametricDistance(pcoords, 1); +} + +template +VTKC_EXEC inline bool cellInside(Line, const CoordType& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + using T = ComponentType; + return component(pcoords, 0) >= T{0} && component(pcoords, 0) <= T{1}; +} + +template +VTKC_EXEC inline vtkc::ErrorCode interpolate( + Line, + const Values& values, + const CoordType& pcoords, + Result&& result) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + using ProcessingType = internal::ClosestFloatType; + using ResultCompType = ComponentType; + + for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c) + { + auto ival = internal::lerp(static_cast(values.getValue(0, c)), + static_cast(values.getValue(1, c)), + static_cast(component(pcoords, 0))); + component(result, c) = static_cast(ival); + } + + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline vtkc::ErrorCode derivative( + Line, + const Points& points, + const Values& values, + const CoordType&, + Result&& dx, + Result&& dy, + Result&& dz) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + using ProcessingType = internal::ClosestFloatType; + using ResultCompType = ComponentType; + + ProcessingType dPt[3] = { + static_cast(points.getValue(1, 0) - points.getValue(0, 0)), + static_cast(points.getValue(1, 1) - points.getValue(0, 1)), + static_cast(points.getValue(1, 2) - points.getValue(0, 2)) }; + + for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c) + { + auto dv = static_cast(values.getValue(1, c) - values.getValue(0, c)); + component(dx, c) = (dPt[0] != 0.0f) ? static_cast(dv/dPt[0]) : ResultCompType{0}; + component(dy, c) = (dPt[1] != 0.0f) ? static_cast(dv/dPt[1]) : ResultCompType{0}; + component(dz, c) = (dPt[2] != 0.0f) ? static_cast(dv/dPt[2]) : ResultCompType{0}; + } + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricToWorld( + Line, + const Points& points, + const PCoordType& pcoords, + WCoordType&& wcoords) noexcept +{ + return interpolate(Line{}, points, pcoords, std::forward(wcoords)); +} + +template +VTKC_EXEC inline vtkc::ErrorCode worldToParametric( + Line, + const Points& points, + const WCoordType& wcoords, + PCoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType); + + using T = ComponentType; + internal::Vector p0(static_cast(points.getValue(0, 0)), + static_cast(points.getValue(0, 1)), + static_cast(points.getValue(0, 2))); + internal::Vector p1(static_cast(points.getValue(1, 0)), + static_cast(points.getValue(1, 1)), + static_cast(points.getValue(1, 2))); + internal::Vector wc(static_cast(component(wcoords, 0)), + static_cast(component(wcoords, 1)), + static_cast(component(wcoords, 2))); + auto v1 = p1 - p0; + auto v2 = wc - p0; + component(pcoords, 0) = internal::dot(v1, v2) / internal::dot(v1, v1); + return ErrorCode::SUCCESS; +} + +} //namespace vtkc + +#endif //vtk_c_Line_h diff --git a/vtkc/Pixel.h b/vtkc/Pixel.h new file mode 100644 index 000000000..8bf42653b --- /dev/null +++ b/vtkc/Pixel.h @@ -0,0 +1,178 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.md 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. +//============================================================================ +#ifndef vtk_c_Pixel_h +#define vtk_c_Pixel_h + +#include +#include + +#include + +namespace vtkc +{ + +class Pixel : public Quad +{ +public: + constexpr VTKC_EXEC Pixel() : Quad(Cell(ShapeId::PIXEL, 4)) {} + constexpr VTKC_EXEC explicit Pixel(const Cell& cell) : Quad(cell) {} +}; + +namespace internal +{ + +template +VTKC_EXEC inline int getPixelSpacing(const Points& points, T spacing[3]) +{ + int zeros = 0; + for (int i = 0; i < 3; ++i) + { + spacing[i] = static_cast(points.getValue(2, i) - points.getValue(0, i)); + if (spacing[i] == T{0}) + { + zeros |= 1 << i; + } + } + return zeros; +} + +} // internal + +template +VTKC_EXEC inline vtkc::ErrorCode derivative( + Pixel, + const Points& points, + const Values& values, + const CoordType& pcoords, + Result&& dx, + Result&& dy, + Result&& dz) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + using ProcessingType = internal::ClosestFloatType; + using ResultCompType = ComponentType; + + ProcessingType spacing[3]; + int zeros = internal::getPixelSpacing(points, spacing); + + for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c) + { + ProcessingType dvdp[2]; + internal::parametricDerivative(Quad{}, values, c, pcoords, dvdp); + switch (zeros) + { + case 1: // yz plane + component(dx, c) = ResultCompType{0}; + component(dy, c) = static_cast(dvdp[0] / spacing[1]); + component(dz, c) = static_cast(dvdp[1] / spacing[2]); + break; + case 2: // xz plane + component(dx, c) = static_cast(dvdp[0] / spacing[0]); + component(dy, c) = ResultCompType{0}; + component(dz, c) = static_cast(dvdp[1] / spacing[2]); + break; + case 4: // xy plane + component(dx, c) = static_cast(dvdp[0] / spacing[0]); + component(dy, c) = static_cast(dvdp[1] / spacing[1]); + component(dz, c) = ResultCompType{0}; + break; + default: + return ErrorCode::DEGENERATE_CELL_DETECTED; + } + } + + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC vtkc::ErrorCode parametricToWorld( + Pixel, + const Points& points, + const PCoordType& pcoords, + WCoordType&& wcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType); + + using T = typename Points::ValueType; + + T spacing[3]; + auto zeros = internal::getPixelSpacing(points, spacing); + + switch (zeros) + { + case 1: // yz plane + component(wcoords, 0) = points.getValue(0, 0); + component(wcoords, 1) = points.getValue(0, 1) + + (spacing[1] * static_cast(component(pcoords, 0))); + component(wcoords, 2) = points.getValue(0, 2) + + (spacing[2] * static_cast(component(pcoords, 1))); + return ErrorCode::SUCCESS; + case 2: // xz plane + component(wcoords, 0) = points.getValue(0, 0) + + (spacing[0] * static_cast(component(pcoords, 0))); + component(wcoords, 1) = points.getValue(0, 1); + component(wcoords, 2) = points.getValue(0, 2) + + (spacing[2] * static_cast(component(pcoords, 1))); + return ErrorCode::SUCCESS; + case 4: // xy plane + component(wcoords, 0) = points.getValue(0, 0) + + (spacing[0] * static_cast(component(pcoords, 0))); + component(wcoords, 1) = points.getValue(0, 1) + + (spacing[1] * static_cast(component(pcoords, 1))); + component(wcoords, 2) = points.getValue(0, 2); + return ErrorCode::SUCCESS; + default: + return ErrorCode::DEGENERATE_CELL_DETECTED; + } +} + +template +VTKC_EXEC vtkc::ErrorCode worldToParametric( + Pixel, + const Points& points, + const WCoordType& wcoords, + PCoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType); + + using T = ComponentType; + + T spacing[3]; + int zeros = internal::getPixelSpacing(points, spacing); + + switch (zeros) + { + case 1: // yz plane + component(pcoords, 0) = static_cast(component(wcoords, 1) - points.getValue(0, 1)) / + spacing[1]; + component(pcoords, 1) = static_cast(component(wcoords, 2) - points.getValue(0, 2)) / + spacing[2]; + return ErrorCode::SUCCESS; + case 2: // xz plane + component(pcoords, 0) = static_cast(component(wcoords, 0) - points.getValue(0, 0)) / + spacing[0]; + component(pcoords, 1) = static_cast(component(wcoords, 2) - points.getValue(0, 2)) / + spacing[2]; + return ErrorCode::SUCCESS; + case 4: // xy plane + component(pcoords, 0) = static_cast(component(wcoords, 0) - points.getValue(0, 0)) / + spacing[0]; + component(pcoords, 1) = static_cast(component(wcoords, 1) - points.getValue(0, 1)) / + spacing[1]; + return ErrorCode::SUCCESS; + default: + return ErrorCode::DEGENERATE_CELL_DETECTED; + } +} + +} //namespace vtkc + +#endif //vtk_c_Pixel_h diff --git a/vtkc/Polygon.h b/vtkc/Polygon.h new file mode 100644 index 000000000..0165c123b --- /dev/null +++ b/vtkc/Polygon.h @@ -0,0 +1,541 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.md 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. +//============================================================================ +#ifndef vtk_c_Polygon_h +#define vtk_c_Polygon_h + +#include +#include +#include +#include + +#include + +namespace vtkc +{ + +/// \c Polygon with 3 and 4 points behave exactly as \c Triangle and \c Quad +/// respectively. For 5 or more points, the points are arranged 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 the 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. +class Polygon : public Cell +{ +public: + constexpr VTKC_EXEC Polygon() : Cell(ShapeId::POLYGON, 3) {} + constexpr VTKC_EXEC explicit Polygon(vtkc::IdComponent numPoints) + : Cell(ShapeId::POLYGON, numPoints) + { + } + constexpr VTKC_EXEC explicit Polygon(const Cell& cell) : Cell(cell) {} +}; + +VTKC_EXEC inline vtkc::ErrorCode validate(Polygon tag) noexcept +{ + if (tag.shape() != ShapeId::POLYGON) + { + return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE; + } + if (tag.numberOfPoints() < 3) + { + return ErrorCode::INVALID_NUMBER_OF_POINTS; + } + + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricCenter(Polygon tag, CoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + switch (tag.numberOfPoints()) + { + case 3: + return parametricCenter(Triangle{}, pcoords); + case 4: + return parametricCenter(Quad{}, pcoords); + default: + component(pcoords, 0) = 0.5f; + component(pcoords, 1) = 0.5f; + return ErrorCode::SUCCESS; + } +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricPoint( + Polygon tag, IdComponent pointId, CoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + if (pointId < 0 || pointId >= tag.numberOfPoints()) + { + return ErrorCode::INVALID_POINT_ID; + } + + switch (tag.numberOfPoints()) + { + case 3: + return parametricPoint(Triangle{}, pointId, pcoords); + case 4: + return parametricPoint(Quad{}, pointId, pcoords); + default: + { + using T = ComponentType; + constexpr double two_pi = 2.0 * 3.14159265359; + auto angle = (static_cast(pointId) * static_cast(two_pi)) / static_cast(tag.numberOfPoints()); + component(pcoords, 0) = 0.5f * (VTKC_MATH_CALL(cos, (angle)) + 1.0f); + component(pcoords, 1) = 0.5f * (VTKC_MATH_CALL(sin, (angle)) + 1.0f); + return ErrorCode::SUCCESS; + } + } +} + +template +VTKC_EXEC inline ComponentType parametricDistance(Polygon tag, const CoordType& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + switch (tag.numberOfPoints()) + { + case 3: + return parametricDistance(Triangle{}, pcoords); + default: + return internal::findParametricDistance(pcoords, 2); + } +} + +template +VTKC_EXEC inline bool cellInside(Polygon tag, const CoordType& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + using T = ComponentType; + switch (tag.numberOfPoints()) + { + case 3: + return cellInside(Triangle{}, pcoords); + case 4: + return cellInside(Quad{}, pcoords); + default: + return (((component(pcoords, 0) - T(0.5f)) * (component(pcoords, 0) - T(0.5f))) + + ((component(pcoords, 1) - T(0.5f)) * (component(pcoords, 1) - T(0.5f)))) <= T(0.25f); + } +} + +namespace internal +{ + +template +VTKC_EXEC inline vtkc::ErrorCode polygonToSubTrianglePCoords( + Polygon tag, + const CoordType& polygonPC, + IdComponent& p0, + IdComponent& p1, + ComponentType trianglePC[2]) noexcept +{ + using T = ComponentType; + + constexpr T epsilon = std::is_same::value ? T(1e-5f) : T(1e-9f); + + // Find the sub-triangle containing pcoords + auto x = component(polygonPC, 1) - T(0.5f); + auto y = component(polygonPC, 0) - T(0.5f); + if (VTKC_MATH_CALL(abs, (x)) < (T(4) * epsilon) && VTKC_MATH_CALL(abs, (y)) < (T(4) * epsilon)) + { + // we are at the center + p0 = 0; + p1 = 1; + trianglePC[0] = trianglePC[1] = T(0); + return ErrorCode::SUCCESS; + } + + constexpr double two_pi = 2.0 * 3.14159265359; + T angle = VTKC_MATH_CALL(atan2, (x), (y)); + if (angle < T(0)) + { + angle += static_cast(two_pi); + } + T deltaAngle = static_cast(two_pi) / static_cast(tag.numberOfPoints()); + + p0 = static_cast(VTKC_MATH_CALL(floor, (angle / deltaAngle))); + p1 = (p0 + 1) % tag.numberOfPoints(); + + // Build triangle with polygon pcoords as its wcoords + T triPts[9] = { T(0), T(0), T(0), T(0), T(0), T(0), T(0), T(0), T(0) }; + VTKC_RETURN_ON_ERROR(parametricCenter(tag, triPts)) + VTKC_RETURN_ON_ERROR(parametricPoint(tag, p0, triPts + 3)) + VTKC_RETURN_ON_ERROR(parametricPoint(tag, p1, triPts + 6)) + + // Find the parametric coord on the triangle + T triWC[3] = { component(polygonPC, 0), component(polygonPC, 1), T(0) }; + VTKC_RETURN_ON_ERROR(worldToParametric(Triangle{}, makeFieldAccessorFlatSOAConst(triPts, 3), triWC, trianglePC)) + + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline typename Values::ValueType polygonInterpolateComponentAtCenter( + Polygon tag, const Values& values, IdComponent comp) noexcept +{ + using T = internal::ClosestFloatType; + + auto weight = T{1} / static_cast(tag.numberOfPoints()); + auto result = static_cast(values.getValue(0, comp)); + for (IdComponent i = 1; i < tag.numberOfPoints(); ++i) + { + result += static_cast(values.getValue(i, comp)); + } + result *= weight; + + return static_cast(result); +} + +} // namespace internal + +template +VTKC_EXEC vtkc::ErrorCode interpolate( + Polygon tag, + const Values& values, + const CoordType& pcoords, + Result&& result) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + switch (tag.numberOfPoints()) + { + case 3: + return interpolate(Triangle{}, values, pcoords, std::forward(result)); + case 4: + return interpolate(Quad{}, values, pcoords, std::forward(result)); + default: + break; + } + + using ResultCompType = ComponentType; + using ProcessingType = internal::ClosestFloatType; + + IdComponent p0, p1; + ComponentType triPc[2]; + VTKC_RETURN_ON_ERROR(internal::polygonToSubTrianglePCoords(tag, pcoords, p0, p1, triPc)) + + // compute polygon interpolation from triangle weights + for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c) + { + ProcessingType triVals[3]; + triVals[0] = static_cast(internal::polygonInterpolateComponentAtCenter(tag, values, c)); + triVals[1] = static_cast(values.getValue(p0, c)); + triVals[2] = static_cast(values.getValue(p1, c)); + ResultCompType val = 0; + VTKC_RETURN_ON_ERROR(interpolate(Triangle{}, makeFieldAccessorNestedSOA(triVals), triPc, &val)) + component(result, c) = val; + } + + return ErrorCode::SUCCESS; +} + +namespace internal +{ +// To find the gradient in a polygon (of 5 or more points), we will extract a small triangle near +// the desired parameteric coordinates (pcoords). We return the field values (outField) and world +// coordinates (outWCoords) for this triangle, which is all that is needed to find the gradient +// in a triangle. +// +// The trangle will be "pointing" away from the center of the polygon, and pcoords will be placed +// at the apex of the triangle. This is because if pcoords is at or near the edge of the polygon, +// we do not want to push any of the points over the edge, and it is not trivial to determine +// exactly where the edge of the polygon is. +template +VTKC_EXEC inline void polygonGetTriangleAroundPCoords( + const CoordType& pcoords, ComponentType pc1[2], ComponentType pc2[2]) noexcept +{ + using T = ComponentType; + + // Find the unit vector pointing from the center of the polygon to pcoords + Vector radialVector(component(pcoords, 0) - 0.5f, component(pcoords, 1) - 0.5f); + auto magSqr = dot(radialVector, radialVector); + if (magSqr > 8.0f * 1e-4f) + { + radialVector /= VTKC_MATH_CALL(sqrt, (magSqr)); + } + else + { + // pcoords is in the center of the polygon. Just point in an arbitrary direction + radialVector[0] = T(1); + radialVector[1] = T(0); + } + + // We want the two points away from pcoords to be back toward the center but moved at 45 degrees + // off the radius. Simple geometry shows us that the (not quite unit) vectors of those two + // directions are (-radialVector[1] - radialVector[0], radialVector[0] - radialVector[1]) and + // (radialVector[1] - radialVector[0], -radialVector[0] - radialVector[1]). + // + // *\ (-radialVector[1], radialVector[0]) // + // | \ // + // | \ (-radialVector[1] - radialVector[0], radialVector[0] - radialVector[1]) // + // | \ // + // +-------* radialVector // + // | / // + // | / (radialVector[1] - radialVector[0], -radialVector[0] - radialVector[1]) // + // | / // + // */ (radialVector[1], -radialVector[0]) // + + // This scaling value is somewhat arbitrary. It is small enough to be "close" to the selected + // point and small enough to be guaranteed to be inside the polygon, but large enough to + // get an accurate gradient. + static constexpr T scale = 0.05f; + + pc1[0] = pcoords[0] + scale * (-radialVector[1] - radialVector[0]); + pc1[1] = pcoords[1] + scale * (radialVector[0] - radialVector[1]); + + pc2[0] = pcoords[0] + scale * (radialVector[1] - radialVector[0]); + pc2[1] = pcoords[1] + scale * (-radialVector[0] - radialVector[1]); +} + +} // namespace internal + +template +VTKC_EXEC inline vtkc::ErrorCode derivative( + Polygon tag, + const Points& points, + const Values& values, + const CoordType& pcoords, + Result&& dx, + Result&& dy, + Result&& dz) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + switch (tag.numberOfPoints()) + { + case 3: + return derivative(Triangle{}, + points, + values, + pcoords, + std::forward(dx), + std::forward(dy), + std::forward(dz)); + case 4: + return derivative(Quad{}, + points, + values, + pcoords, + std::forward(dx), + std::forward(dy), + std::forward(dz)); + default: + break; + } + + using ResultCompType = ComponentType; + using ProcessingType = internal::ClosestFloatType; + + // Get the parametric coordinates of a small triangle, with pcoords as one of the vertices + ComponentType ptPc1[2], ptPc2[2]; + internal::polygonGetTriangleAroundPCoords(pcoords, ptPc1, ptPc2); + + // Compute world coordinates of the points of the triangle + internal::Vector triPts[3]; + VTKC_RETURN_ON_ERROR(interpolate(tag, points, pcoords, triPts[0])) + VTKC_RETURN_ON_ERROR(interpolate(tag, points, ptPc1, triPts[1])) + VTKC_RETURN_ON_ERROR(interpolate(tag, points, ptPc2, triPts[2])) + + // Compute the derivative on the triangle + //---------------------------------------- + // 2-D coordinate system on the triangle's plane + internal::Space2D triSpace(triPts[0], triPts[1], triPts[2]); + internal::Vector pts2d[3]; + for (int i = 0; i < 3; ++i) + { + pts2d[i] = triSpace.to2DPoint(triPts[i]); + } + + // pre-compute once + internal::Matrix jacobian; + internal::jacobian2D(Triangle{}, makeFieldAccessorNestedSOA(pts2d, 2), nullptr, jacobian); + internal::Matrix invJacobian; + VTKC_RETURN_ON_ERROR(internal::matrixInverse(jacobian, invJacobian)) + + // Compute sub-triangle information of the three vertices of the derivation triangle to + // reduce the amount of redundant computations in the loop. + IdComponent subP1P2[3][2]; + ComponentType pcs[3][2]; + internal::polygonToSubTrianglePCoords(tag, pcoords, subP1P2[0][0], subP1P2[0][1], pcs[0]); + internal::polygonToSubTrianglePCoords(tag, ptPc1, subP1P2[1][0], subP1P2[1][1], pcs[1]); + internal::polygonToSubTrianglePCoords(tag, ptPc2, subP1P2[2][0], subP1P2[2][1], pcs[2]); + + for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c) + { + // Interpolate component values at the vertices of the derivation triangle. + auto vCenter = static_cast(internal::polygonInterpolateComponentAtCenter(tag, values, c)); + ProcessingType triVals[3]; + for (int i = 0; i < 3; ++i) + { + ProcessingType field[3] = {vCenter, + static_cast(values.getValue(subP1P2[i][0], c)), + static_cast(values.getValue(subP1P2[i][1], c))}; + VTKC_RETURN_ON_ERROR(interpolate(Triangle{}, makeFieldAccessorNestedSOA(field), pcs[i], triVals + i)) + } + + // Compute derivative in the triangle + internal::Vector dvdp; + parametricDerivative(Triangle{}, makeFieldAccessorNestedSOA(triVals), 0, nullptr, dvdp); + auto d2D = matrixMultiply(dvdp, invJacobian); + auto d3D = triSpace.to3DVec(d2D); + + component(dx, c) = static_cast(d3D[0]); + component(dy, c) = static_cast(d3D[1]); + component(dz, c) = static_cast(d3D[2]); + } + + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricToWorld( + Polygon tag, + const Points& points, + const PCoordType& pcoords, + WCoordType&& wcoords) noexcept +{ + return interpolate(tag, points, pcoords, std::forward(wcoords)); +} + +template +VTKC_EXEC inline vtkc::ErrorCode worldToParametric( + Polygon tag, + const Points& points, + const WCoordType& wcoords, + PCoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType); + + switch (tag.numberOfPoints()) + { + case 3: + return worldToParametric(Triangle{}, points, wcoords, pcoords); + case 4: + return worldToParametric(Quad{}, points, wcoords, pcoords); + default: + break; + } + + using T = ComponentType; + auto numPoints = tag.numberOfPoints(); + + // Find the position of the center point. + internal::Vector wcoordCenter{T(0), T(0), T(0)}; + for (IdComponent pointIndex = 0; pointIndex < numPoints; ++pointIndex) + { + wcoordCenter[0] += points.getValue(pointIndex, 0); + wcoordCenter[1] += points.getValue(pointIndex, 1); + wcoordCenter[2] += points.getValue(pointIndex, 2); + } + wcoordCenter /= static_cast(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. + internal::Vector polygonNormal; + { + internal::Vector v1p1, v1p2; + points.getTuple(0, v1p1); + points.getTuple(numPoints / 3, v1p2); + + internal::Vector v2p1, v2p2; + points.getTuple(1, v2p1); + points.getTuple(2 * numPoints / 3, v2p2); + + polygonNormal = internal::cross(v1p2 - v1p1, v2p2 - v2p1); + } + + // 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. + internal::Vector wc{ component(wcoords, 0), component(wcoords, 1), component(wcoords, 2) }; + IdComponent firstPointIndex; + IdComponent secondPointIndex = 0; + internal::Vector firstPoint, secondPoint; + bool foundTriangle = false; + for (firstPointIndex = 0; firstPointIndex < numPoints - 1; ++firstPointIndex) + { + points.getTuple(firstPointIndex, firstPoint); + auto vecInPlane = firstPoint - wcoordCenter; + + auto planeNormal = internal::cross(polygonNormal, vecInPlane); + auto planeOffset = internal::dot(planeNormal, wcoordCenter); + if (internal::dot(planeNormal, wc) < planeOffset) + { + // wcoords on wrong side of plane, thus outside of triangle + continue; + } + + secondPointIndex = firstPointIndex + 1; + points.getTuple(secondPointIndex, secondPoint); + vecInPlane = secondPoint - wcoordCenter; + + planeNormal = internal::cross(polygonNormal, vecInPlane); + planeOffset = internal::dot(planeNormal, wcoordCenter); + if (internal::dot(planeNormal, wc) > 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; + points.getTuple(firstPointIndex, firstPoint); + secondPointIndex = 0; + points.getTuple(secondPointIndex, secondPoint); + } + + // 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. + internal::Vector triangleWCoords[3] = { wcoordCenter, firstPoint, secondPoint }; + internal::Vector trianglePCoords; + VTKC_RETURN_ON_ERROR(worldToParametric( + Triangle{}, makeFieldAccessorNestedSOA(triangleWCoords, 3), wc, trianglePCoords)) + + // 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 parametricToWorld by using the + // polygon parametric coordinates as a proxy for world coordinates. + VTKC_RETURN_ON_ERROR(parametricCenter(tag, triangleWCoords[0])) + VTKC_RETURN_ON_ERROR(parametricPoint(tag, firstPointIndex, triangleWCoords[1])) + VTKC_RETURN_ON_ERROR(parametricPoint(tag, secondPointIndex, triangleWCoords[2])) + triangleWCoords[0][2] = triangleWCoords[1][2] = triangleWCoords[2][2] = T(0); + VTKC_RETURN_ON_ERROR( + parametricToWorld(Triangle{}, makeFieldAccessorNestedSOA(triangleWCoords, 3), trianglePCoords, wc)) + + component(pcoords, 0) = static_cast>(wc[0]); + component(pcoords, 1) = static_cast>(wc[1]); + + return ErrorCode::SUCCESS; +} + +} //namespace vtkc + +#endif //vtk_c_Polygon_h diff --git a/vtkc/Pyramid.h b/vtkc/Pyramid.h new file mode 100644 index 000000000..80f72ba92 --- /dev/null +++ b/vtkc/Pyramid.h @@ -0,0 +1,293 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.md 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. +//============================================================================ +#ifndef vtk_c_Pyramid_h +#define vtk_c_Pyramid_h + +#include +#include + +#include + +namespace vtkc +{ + +class Pyramid : public Cell +{ +public: + constexpr VTKC_EXEC Pyramid() : Cell(ShapeId::PYRAMID, 5) {} + constexpr VTKC_EXEC explicit Pyramid(const Cell& cell) : Cell(cell) {} +}; + +VTKC_EXEC inline vtkc::ErrorCode validate(Pyramid tag) noexcept +{ + if (tag.shape() != ShapeId::PYRAMID) + { + return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE; + } + if (tag.numberOfPoints() != 5) + { + return ErrorCode::INVALID_NUMBER_OF_POINTS; + } + + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricCenter(Pyramid, CoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + component(pcoords, 0) = 0.5f; + component(pcoords, 1) = 0.5f; + component(pcoords, 2) = 0.2f; + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricPoint( + Pyramid, IdComponent pointId, CoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + switch (pointId) + { + case 0: + component(pcoords, 0) = 0.0f; + component(pcoords, 1) = 0.0f; + component(pcoords, 2) = 0.0f; + break; + case 1: + component(pcoords, 0) = 1.0f; + component(pcoords, 1) = 0.0f; + component(pcoords, 2) = 0.0f; + break; + case 2: + component(pcoords, 0) = 1.0f; + component(pcoords, 1) = 1.0f; + component(pcoords, 2) = 0.0f; + break; + case 3: + component(pcoords, 0) = 0.0f; + component(pcoords, 1) = 1.0f; + component(pcoords, 2) = 0.0f; + break; + case 4: + component(pcoords, 0) = 0.5f; + component(pcoords, 1) = 0.5f; + component(pcoords, 2) = 1.0f; + break; + default: + return ErrorCode::INVALID_POINT_ID; + } + + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline ComponentType parametricDistance(Pyramid, const CoordType& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + return internal::findParametricDistance(pcoords, 3); +} + +template +VTKC_EXEC inline bool cellInside(Pyramid, const CoordType& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + using T = ComponentType; + return component(pcoords, 0) >= T{0} && component(pcoords, 0) <= T{1} && + component(pcoords, 1) >= T{0} && component(pcoords, 1) <= T{1} && + component(pcoords, 2) >= T{0} && component(pcoords, 2) <= T{1}; +} + +template +VTKC_EXEC inline vtkc::ErrorCode interpolate( + Pyramid, + const Values& values, + const CoordType& pcoords, + Result&& result) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + using T = internal::ClosestFloatType; + + for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c) + { + auto baseV0 = internal::lerp(static_cast(values.getValue(0, c)), + static_cast(values.getValue(1, c)), + static_cast(component(pcoords, 0))); + auto baseV1 = internal::lerp(static_cast(values.getValue(3, c)), + static_cast(values.getValue(2, c)), + static_cast(component(pcoords, 0))); + + auto baseV = internal::lerp(baseV0, baseV1, static_cast(component(pcoords, 1))); + + auto v = internal::lerp(baseV, + static_cast(values.getValue(4, c)), + static_cast(component(pcoords, 2))); + + component(result, c) = static_cast>(v); + } + + return ErrorCode::SUCCESS; +} + +namespace internal +{ + +template +VTKC_EXEC inline void parametricDerivative(Pyramid, + const Values& values, + IdComponent comp, + const CoordType& pcoords, + Result&& result) noexcept +{ + using T = internal::ClosestFloatType; + T p0 = static_cast(component(pcoords, 0)); + T p1 = static_cast(component(pcoords, 1)); + T p2 = static_cast(component(pcoords, 2)); + T rm = T{1} - p0; + T sm = T{1} - p1; + T tm = T{1} - p2; + + T dr = (static_cast(values.getValue(0, comp)) * -sm * tm) + + (static_cast(values.getValue(1, comp)) * sm * tm) + + (static_cast(values.getValue(2, comp)) * p1 * tm) + + (static_cast(values.getValue(3, comp)) * -p1 * tm); + + T ds = (static_cast(values.getValue(0, comp)) * -rm * tm) + + (static_cast(values.getValue(1, comp)) * -p0 * tm) + + (static_cast(values.getValue(2, comp)) * p0 * tm) + + (static_cast(values.getValue(3, comp)) * rm * tm); + + T dt = (static_cast(values.getValue(0, comp)) * -rm * sm) + + (static_cast(values.getValue(1, comp)) * -p0 * sm) + + (static_cast(values.getValue(2, comp)) * -p0 * p1) + + (static_cast(values.getValue(3, comp)) * -rm * p1) + + (static_cast(values.getValue(4, comp))); + + component(result, 0) = static_cast>(dr); + component(result, 1) = static_cast>(ds); + component(result, 2) = static_cast>(dt); +} + +} // internal + +template +VTKC_EXEC inline vtkc::ErrorCode derivative( + Pyramid, + const Points& points, + const Values& values, + const CoordType& pcoords, + Result&& dx, + Result&& dy, + Result&& dz) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + using ProcessingType = internal::ClosestFloatType; + using ResultCompType = ComponentType; + + if (component(pcoords, 2) > ComponentType(.999f)) + { + // If we are at the apex of the pyramid we need to do something special. + // As we approach the apex, the derivatives of the parametric shape + // functions in x and y go to 0 while the inverse of the Jacobian + // also goes to 0. This results in 0/0 but using l'Hopital's rule + // we could actually compute the value of the limit, if we had a + // functional expression to compute the gradient. We're on a computer + // so we don't but we can cheat and do a linear extrapolation of the + // derivatives which really ends up as the same thing. + internal::Matrix j, ij1, ij2; + + const ComponentType pc1[3] = {0.5f, 0.5f, (2.0f * 0.998f) - component(pcoords, 2)}; + internal::jacobian3D(Pyramid{}, points, pc1, j); + VTKC_RETURN_ON_ERROR(internal::matrixInverse(j, ij1)) + + const ComponentType pc2[3] = {0.5f, 0.5f, 0.998f}; + internal::jacobian3D(Pyramid{}, points, pc2, j); + VTKC_RETURN_ON_ERROR(internal::matrixInverse(j, ij2)) + + for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c) + { + internal::Vector dvdp; + + internal::parametricDerivative(Pyramid{}, values, c, pc1, dvdp); + auto d1 = matrixMultiply(dvdp, ij1); + + internal::parametricDerivative(Pyramid{}, values, c, pc2, dvdp); + auto d2 = matrixMultiply(dvdp, ij2); + + component(dx, c) = static_cast((d2[0] * 2.0f) - d1[0]); + component(dy, c) = static_cast((d2[1] * 2.0f) - d1[1]); + component(dz, c) = static_cast((d2[2] * 2.0f) - d1[2]); + } + } + else + { + return internal::derivative3D(Pyramid{}, + points, + values, + pcoords, + std::forward(dx), + std::forward(dy), + std::forward(dz)); + } + + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricToWorld( + Pyramid, + const Points& points, + const PCoordType& pcoords, + WCoordType&& wcoords) noexcept +{ + return interpolate(Pyramid{}, points, pcoords, std::forward(wcoords)); +} + +template +VTKC_EXEC inline vtkc::ErrorCode worldToParametric( + Pyramid, + const Points& points, + const WCoordType& wcoords, + PCoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType); + + using TIn = typename Points::ValueType; + using TOut = ComponentType; + + internal::Vector wcVec{component(wcoords, 0), component(wcoords, 1), component(wcoords, 2)}; + + // Newton's method fails if the wcoord is too close to the apex. Just return the pcoords at the + // apex for those cases. + internal::Vector pcBaseCenter(0.5f, 0.5f, 0.0f); + internal::Vector apex, wcBaseCenter; + points.getTuple(4, apex); + VTKC_RETURN_ON_ERROR(parametricToWorld(Pyramid{}, points, pcBaseCenter, wcBaseCenter)) + auto apexToBase = wcBaseCenter - apex; + auto apexToWc = wcVec - apex; + auto dist2ApexToBase = internal::dot(apexToBase, apexToBase); + auto dist2ApexToWC = internal::dot(apexToWc, apexToWc); + if (dist2ApexToWC <= (1e-6f * dist2ApexToBase)) + { + return parametricPoint(Pyramid{}, 4, pcoords); + } + + return internal::worldToParametric3D( + Pyramid{}, points, wcoords, std::forward(pcoords)); +} + +} // vtkc + +#endif // vtk_c_Pyramid_h diff --git a/vtkc/Quad.h b/vtkc/Quad.h new file mode 100644 index 000000000..861d28ed3 --- /dev/null +++ b/vtkc/Quad.h @@ -0,0 +1,195 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.md 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. +//============================================================================ +#ifndef vtk_c_Quad_h +#define vtk_c_Quad_h + +#include +#include + +#include + +namespace vtkc +{ + +class Quad : public Cell +{ +public: + constexpr VTKC_EXEC Quad() : Cell(ShapeId::QUAD, 4) {} + constexpr VTKC_EXEC explicit Quad(const Cell& cell) : Cell(cell) {} +}; + +VTKC_EXEC inline vtkc::ErrorCode validate(Quad tag) noexcept +{ + if (tag.shape() != ShapeId::QUAD && tag.shape() != ShapeId::PIXEL) + { + return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE; + } + if (tag.numberOfPoints() != 4) + { + return ErrorCode::INVALID_NUMBER_OF_POINTS; + } + + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricCenter(Quad, CoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + component(pcoords, 0) = 0.5f; + component(pcoords, 1) = 0.5f; + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricPoint( + Quad, IdComponent pointId, CoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + switch (pointId) + { + case 0: + component(pcoords, 0) = 0.0f; + component(pcoords, 1) = 0.0f; + break; + case 1: + component(pcoords, 0) = 1.0f; + component(pcoords, 1) = 0.0f; + break; + case 2: + component(pcoords, 0) = 1.0f; + component(pcoords, 1) = 1.0f; + break; + case 3: + component(pcoords, 0) = 0.0f; + component(pcoords, 1) = 1.0f; + break; + default: + return ErrorCode::INVALID_POINT_ID; + } + + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline ComponentType parametricDistance(Quad, const CoordType& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + return internal::findParametricDistance(pcoords, 2); +} + +template +VTKC_EXEC inline bool cellInside(Quad, const CoordType& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + using T = ComponentType; + return component(pcoords, 0) >= T{0} && component(pcoords, 0) <= T{1} && + component(pcoords, 1) >= T{0} && component(pcoords, 1) <= T{1}; +} + +template +VTKC_EXEC inline vtkc::ErrorCode interpolate( + Quad, + const Values& values, + const CoordType& pcoords, + Result&& result) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + using T = internal::ClosestFloatType; + + for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c) + { + auto v0 = internal::lerp(static_cast(values.getValue(0, c)), + static_cast(values.getValue(1, c)), + static_cast(component(pcoords, 0))); + auto v1 = internal::lerp(static_cast(values.getValue(3, c)), + static_cast(values.getValue(2, c)), + static_cast(component(pcoords, 0))); + auto v = internal::lerp(v0, v1, static_cast(component(pcoords, 1))); + component(result, c) = static_cast>(v); + } + + return ErrorCode::SUCCESS; +} + +namespace internal +{ + +template +VTKC_EXEC inline void parametricDerivative( + Quad, const Values& values, IdComponent comp, const CoordType& pcoords, Result&& result) noexcept +{ + using T = internal::ClosestFloatType; + T p0 = static_cast(component(pcoords, 0)); + T p1 = static_cast(component(pcoords, 1)); + T rm = T{1} - p0; + T sm = T{1} - p1; + + T dr = (static_cast(values.getValue(0, comp)) * -sm) + + (static_cast(values.getValue(1, comp)) * sm) + + (static_cast(values.getValue(2, comp)) * p1) + + (static_cast(values.getValue(3, comp)) * -p1); + T ds = (static_cast(values.getValue(0, comp)) * -rm) + + (static_cast(values.getValue(1, comp)) * -p0) + + (static_cast(values.getValue(2, comp)) * p0) + + (static_cast(values.getValue(3, comp)) * rm); + + component(result, 0) = static_cast>(dr); + component(result, 1) = static_cast>(ds); +} + +} // internal + +template +VTKC_EXEC inline vtkc::ErrorCode derivative( + Quad, + const Points& points, + const Values& values, + const CoordType& pcoords, + Result&& dx, + Result&& dy, + Result&& dz) noexcept +{ + return internal::derivative2D(Quad{}, + points, + values, + pcoords, + std::forward(dx), + std::forward(dy), + std::forward(dz)); +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricToWorld( + Quad, + const Points& points, + const PCoordType& pcoords, + WCoordType&& wcoords) noexcept +{ + return interpolate(Quad{}, points, pcoords, std::forward(wcoords)); +} + +template +VTKC_EXEC inline vtkc::ErrorCode worldToParametric( + Quad, + const Points& points, + const WCoordType& wcoords, + PCoordType&& pcoords) noexcept +{ + return internal::worldToParametric2D(Quad{}, points, wcoords, std::forward(pcoords)); +} + +} //namespace vtkc + +#endif //vtk_c_Quad_h diff --git a/vtkc/Shapes.h b/vtkc/Shapes.h new file mode 100644 index 000000000..c56476e6b --- /dev/null +++ b/vtkc/Shapes.h @@ -0,0 +1,142 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.md 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. +//============================================================================ +#ifndef vtk_c_Shapes_h +#define vtk_c_Shapes_h + +#include + +#include + +namespace vtkc +{ + +enum ShapeId : IdShape +{ + // Linear cells + EMPTY = 0, + VERTEX = 1, + //POLY_VERTEX = 2, + LINE = 3, + //POLY_LINE = 4, + TRIANGLE = 5, + //TRIANGLE_STRIP = 6, + POLYGON = 7, + PIXEL = 8, + QUAD = 9, + TETRA = 10, + VOXEL = 11, + HEXAHEDRON = 12, + WEDGE = 13, + PYRAMID = 14, + + NUMBER_OF_CELL_SHAPES +}; + +class Cell +{ +public: + constexpr VTKC_EXEC Cell() : Shape(ShapeId::EMPTY), NumberOfPoints(0) {} + constexpr VTKC_EXEC Cell(IdShape shape, IdComponent numberOfPoints) + : Shape(shape), NumberOfPoints(numberOfPoints) + { + } + + constexpr VTKC_EXEC IdShape shape() const noexcept { return this->Shape; } + constexpr VTKC_EXEC IdComponent numberOfPoints() const noexcept { return this->NumberOfPoints; } + +protected: + IdShape Shape; + IdComponent NumberOfPoints; +}; + +/// \brief Check if a shape id is valid +/// \param[in] shapeId The id to check. +/// \return true if id is a valid shape id. +/// +constexpr inline VTKC_EXEC bool isValidShape(IdShape shapeId) +{ + return (shapeId >= ShapeId::EMPTY) && (shapeId < ShapeId::NUMBER_OF_CELL_SHAPES); +} + +/// \brief Returns the dimensionality of a cell +/// \param[in] shapeId The shape id of the cell. +/// \return The dimensionality of the cell, -1 for invalid shapes. +/// +inline VTKC_EXEC int dimension(IdShape shapeId) +{ + switch (shapeId) + { + case VERTEX: + return 0; + case LINE: + return 1; + case TRIANGLE: + case POLYGON: + case PIXEL: + case QUAD: + return 2; + case TETRA: + case VOXEL: + case HEXAHEDRON: + case WEDGE: + case PYRAMID: + return 3; + case EMPTY: + default: + return -1; + } +} + +/// \brief Returns the dimensionality of a cell +/// \param[in] cell The cell. +/// \return The dimensionality of the cell, -1 for invalid shapes. +/// +inline VTKC_EXEC int dimension(Cell cell) +{ + return dimension(cell.shape()); +} + +// forward declare cell tags +class Vertex; +class Line; +class Triangle; +class Polygon; +class Pixel; +class Quad; +class Tetra; +class Voxel; +class Hexahedron; +class Wedge; +class Pyramid; + +} //namespace vtkc + +#define vtkcGenericCellShapeMacroCase(cellId, cell, call) \ + case cellId: \ + { \ + using CellTag = cell; \ + call; \ + } \ + break + +#define vtkcGenericCellShapeMacro(call) \ + vtkcGenericCellShapeMacroCase(vtkc::ShapeId::VERTEX, vtkc::Vertex, call); \ + vtkcGenericCellShapeMacroCase(vtkc::ShapeId::LINE, vtkc::Line, call); \ + vtkcGenericCellShapeMacroCase(vtkc::ShapeId::TRIANGLE, vtkc::Triangle, call); \ + vtkcGenericCellShapeMacroCase(vtkc::ShapeId::POLYGON, vtkc::Polygon, call); \ + vtkcGenericCellShapeMacroCase(vtkc::ShapeId::PIXEL, vtkc::Pixel, call); \ + vtkcGenericCellShapeMacroCase(vtkc::ShapeId::QUAD, vtkc::Quad, call); \ + vtkcGenericCellShapeMacroCase(vtkc::ShapeId::TETRA, vtkc::Tetra, call); \ + vtkcGenericCellShapeMacroCase(vtkc::ShapeId::VOXEL, vtkc::Voxel, call); \ + vtkcGenericCellShapeMacroCase(vtkc::ShapeId::HEXAHEDRON, vtkc::Hexahedron, call); \ + vtkcGenericCellShapeMacroCase(vtkc::ShapeId::WEDGE, vtkc::Wedge, call); \ + vtkcGenericCellShapeMacroCase(vtkc::ShapeId::PYRAMID, vtkc::Pyramid, call) + +#endif //vtk_c_Shapes_h diff --git a/vtkc/Tetra.h b/vtkc/Tetra.h new file mode 100644 index 000000000..5bb326302 --- /dev/null +++ b/vtkc/Tetra.h @@ -0,0 +1,208 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.md 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. +//============================================================================ +#ifndef vtk_c_Tetra_h +#define vtk_c_Tetra_h + +#include +#include + +#include + +namespace vtkc +{ + +class Tetra : public Cell +{ +public: + constexpr VTKC_EXEC Tetra() : Cell(ShapeId::TETRA, 4) {} + constexpr VTKC_EXEC explicit Tetra(const Cell& cell) : Cell(cell) {} +}; + +VTKC_EXEC inline vtkc::ErrorCode validate(Tetra tag) noexcept +{ + if (tag.shape() != ShapeId::TETRA) + { + return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE; + } + if (tag.numberOfPoints() != 4) + { + return ErrorCode::INVALID_NUMBER_OF_POINTS; + } + + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricCenter(Tetra, CoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + component(pcoords, 0) = 0.25f; + component(pcoords, 1) = 0.25f; + component(pcoords, 2) = 0.25f; + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricPoint( + Tetra, IdComponent pointId, CoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + switch (pointId) + { + case 0: + component(pcoords, 0) = 0.0f; + component(pcoords, 1) = 0.0f; + component(pcoords, 2) = 0.0f; + break; + case 1: + component(pcoords, 0) = 1.0f; + component(pcoords, 1) = 0.0f; + component(pcoords, 2) = 0.0f; + break; + case 2: + component(pcoords, 0) = 0.0f; + component(pcoords, 1) = 1.0f; + component(pcoords, 2) = 0.0f; + break; + case 3: + component(pcoords, 0) = 0.0f; + component(pcoords, 1) = 0.0f; + component(pcoords, 2) = 1.0f; + break; + default: + return ErrorCode::INVALID_POINT_ID; + } + + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline ComponentType parametricDistance(Tetra, const CoordType& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + ComponentType weights[4]; + weights[0] = ComponentType{1} - component(pcoords, 0) - component(pcoords, 1) - component(pcoords, 2); + weights[1] = component(pcoords, 0); + weights[2] = component(pcoords, 1); + weights[3] = component(pcoords, 2); + interpolationWeights(Tetra{}, pcoords, weights); + return internal::findParametricDistance(weights, 4); +} + +template +VTKC_EXEC inline bool cellInside(Tetra, const CoordType& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + using T = ComponentType; + return component(pcoords, 0) >= T{0} && + component(pcoords, 1) >= T{0} && + component(pcoords, 2) >= T{0} && + (component(pcoords, 0) + component(pcoords, 1) + component(pcoords, 2)) <= T{1}; +} + +template +VTKC_EXEC inline void interpolationWeights(Tetra, const CoordType& pcoords, T weights[4]) noexcept +{ + weights[1] = static_cast(component(pcoords, 0)); + weights[2] = static_cast(component(pcoords, 1)); + weights[3] = static_cast(component(pcoords, 2)); + weights[0] = T(1) - weights[1] - weights[2] - weights[3]; +} + +template +VTKC_EXEC inline vtkc::ErrorCode interpolate( + Tetra, + const Values& values, + const CoordType& pcoords, + Result&& result) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + using T = internal::ClosestFloatType; + + auto w0 = T(1) - static_cast(component(pcoords, 0) + component(pcoords, 1) + component(pcoords, 2)); + auto w1 = static_cast(component(pcoords, 0)); + auto w2 = static_cast(component(pcoords, 1)); + auto w3 = static_cast(component(pcoords, 2)); + + for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c) + { + auto v = static_cast(values.getValue(0, c)) * w0 + + static_cast(values.getValue(1, c)) * w1 + + static_cast(values.getValue(2, c)) * w2 + + static_cast(values.getValue(3, c)) * w3; + component(result, c) = static_cast>(v); + } + return ErrorCode::SUCCESS; +} + +namespace internal +{ + +template +VTKC_EXEC inline void parametricDerivative( + Tetra, const Values& values, IdComponent comp, const CoordType&, Result&& result) noexcept +{ + component(result, 0) = static_cast>(values.getValue(1, comp) - + values.getValue(0, comp)); + component(result, 1) = static_cast>(values.getValue(2, comp) - + values.getValue(0, comp)); + component(result, 2) = static_cast>(values.getValue(3, comp) - + values.getValue(0, comp)); +} + +} // internal + +template +VTKC_EXEC inline vtkc::ErrorCode derivative( + Tetra, + const Points& points, + const Values& values, + const CoordType& pcoords, + Result&& dx, + Result&& dy, + Result&& dz) noexcept +{ + return internal::derivative3D(Tetra{}, + points, + values, + pcoords, + std::forward(dx), + std::forward(dy), + std::forward(dz)); +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricToWorld( + Tetra, + const Points& points, + const PCoordType& pcoords, + WCoordType&& wcoords) noexcept +{ + return interpolate(Tetra{}, points, pcoords, std::forward(wcoords)); +} + +template +VTKC_EXEC inline vtkc::ErrorCode worldToParametric( + Tetra, + const Points& points, + const WCoordType& wcoords, + PCoordType&& pcoords) noexcept +{ + return internal::worldToParametric3D(Tetra{}, points, wcoords, std::forward(pcoords)); +} + +} // vtkc + +#endif // vtk_c_Tetra_h diff --git a/vtkc/Triangle.h b/vtkc/Triangle.h new file mode 100644 index 000000000..e05f5c4c8 --- /dev/null +++ b/vtkc/Triangle.h @@ -0,0 +1,271 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.md 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. +//============================================================================ +#ifndef vtk_c_Triangle_h +#define vtk_c_Triangle_h + +#include +#include + +#include + +namespace vtkc +{ + +class Triangle : public Cell +{ +public: + constexpr VTKC_EXEC Triangle() : Cell(ShapeId::TRIANGLE, 3) {} + constexpr VTKC_EXEC explicit Triangle(const Cell& cell) : Cell(cell) {} +}; + +VTKC_EXEC inline vtkc::ErrorCode validate(Triangle tag) noexcept +{ + if (tag.shape() != ShapeId::TRIANGLE) + { + return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE; + } + if (tag.numberOfPoints() != 3) + { + return ErrorCode::INVALID_NUMBER_OF_POINTS; + } + + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricCenter(Triangle, CoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + using T = ComponentType; + component(pcoords, 0) = T(1)/T(3); + component(pcoords, 1) = T(1)/T(3); + component(pcoords, 2) = T(0); + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricPoint( + Triangle, IdComponent pointId, CoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + component(pcoords, 2) = 0.0f; + switch (pointId) + { + case 0: + component(pcoords, 0) = 0.0f; + component(pcoords, 1) = 0.0f; + break; + case 1: + component(pcoords, 0) = 1.0f; + component(pcoords, 1) = 0.0f; + break; + case 2: + component(pcoords, 0) = 0.0f; + component(pcoords, 1) = 1.0f; + break; + default: + return ErrorCode::INVALID_POINT_ID; + } + + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline ComponentType parametricDistance(Triangle, const CoordType& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + ComponentType weights[3]; + weights[0] = ComponentType{1} - component(pcoords, 0) - component(pcoords, 1); + weights[1] = component(pcoords, 0); + weights[2] = component(pcoords, 1); + return internal::findParametricDistance(weights, 3); +} + +template +VTKC_EXEC inline bool cellInside(Triangle, const CoordType& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + using T = ComponentType; + return component(pcoords, 0) >= T{0} && component(pcoords, 1) >= T{0} && + (component(pcoords, 0) + component(pcoords, 1)) <= T{1}; +} + +template +VTKC_EXEC inline vtkc::ErrorCode interpolate( + Triangle, + const Values& values, + const CoordType& pcoords, + Result&& result) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + using T = internal::ClosestFloatType; + + auto w0 = T(1) - static_cast(component(pcoords, 0) + component(pcoords, 1)); + auto w1 = static_cast(component(pcoords, 0)); + auto w2 = static_cast(component(pcoords, 1)); + + for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c) + { + auto v = static_cast(values.getValue(0, c)) * w0 + + static_cast(values.getValue(1, c)) * w1 + + static_cast(values.getValue(2, c)) * w2; + component(result, c) = static_cast>(v); + } + return ErrorCode::SUCCESS; +} + +namespace internal +{ + +template +VTKC_EXEC inline void parametricDerivative( + Triangle, const Values& values, IdComponent comp, const CoordType&, Result&& result) noexcept +{ + component(result, 0) = static_cast>(values.getValue(1, comp) - + values.getValue(0, comp)); + component(result, 1) = static_cast>(values.getValue(2, comp) - + values.getValue(0, comp)); +} + +} // internal + +template +VTKC_EXEC inline vtkc::ErrorCode derivative( + Triangle, + const Points& points, + const Values& values, + const CoordType& pcoords, + Result&& dx, + Result&& dy, + Result&& dz) noexcept +{ + return internal::derivative2D(Triangle{}, + points, + values, + pcoords, + std::forward(dx), + std::forward(dy), + std::forward(dz)); +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricToWorld( + Triangle, + const Points& points, + const PCoordType& pcoords, + WCoordType&& wcoords) noexcept +{ + return interpolate(Triangle{}, points, pcoords, std::forward(wcoords)); +} + +//----------------------------------------------------------------------------- +// The following implementation is lifted from VTK-m +//----------------------------------------------------------------------------- +// We will solve the world to parametric coordinates problem geometrically. +// Consider the parallelogram formed by wcoords and p0 of the triangle and +// the two adjacent edges. This parallelogram is equivalent to the +// axis-aligned rectangle anchored at the origin of parametric space. +// +// p2 |\ (1,0) |\ // +// | \ | \ // +// | \ | \ // +// | \ | \ // +// | \ | \ // +// | \ | (u,v) \ // +// | --- \ |-------* \ // +// | ---*wcoords | | \ // +// | | \ | | \ // +// p0 *--- | \ (0,0) *------------------\ (1,0) // +// ---| \ // +// x-- \ // +// --- \ // +// ---\ p1 // +// +// In this diagram, the distance between p0 and the point marked x divided by +// the length of the edge it is on is equal, by proportionality, to the u +// parametric coordinate. (The v coordinate follows the other edge +// accordingly.) Thus, if we can find the intersection at x (or more +// specifically the distance between p0 and x), then we can find that +// parametric coordinate. +// +// Because the triangle is in 3-space, we are actually going to intersect the +// edge with a plane that is parallel to the opposite edge of p0 and +// perpendicular to the triangle. This is partially because it is easy to +// find the intersection between a plane and a line and partially because the +// computation will work for points not on the plane. (The result is +// equivalent to a point projected on the plane.) +// +// First, we define an implicit plane as: +// +// Dot((p - wcoords), planeNormal) = 0 +// +// where planeNormal is the normal to the plane (easily computed from the +// triangle), and p is any point in the plane. Next, we define the parametric +// form of the line: +// +// p(d) = (p1 - p0)d + p0 +// +// Where d is the fraction of distance from p0 toward p1. Note that d is +// actually equal to the parametric coordinate we are trying to find. Once we +// compute it, we are done. We can skip the part about finding the actual +// coordinates of the intersection. +// +// Solving for the intersection is as simple as substituting the line's +// definition of p(d) into p for the plane equation. With some basic algebra +// you get: +// +// d = Dot((wcoords - p0), planeNormal)/Dot((p1-p0), planeNormal) +// +// From here, the u coordinate is simply d. The v coordinate follows +// similarly. +// +template +VTKC_EXEC inline vtkc::ErrorCode worldToParametric( + Triangle, + const Points& points, + const WCoordType& wcoords, + PCoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType); + + using TIn = typename Points::ValueType; + using TOut = ComponentType; + + internal::Vector pts[3]; + for (int i = 0; i < 3; ++i) + { + points.getTuple(i, pts[i]); + } + + internal::Vector wc(component(wcoords, 0), component(wcoords, 1), component(wcoords, 2)); + + auto triangleNormal = internal::cross(pts[1] - pts[0], pts[2] - pts[0]); + for (IdComponent i = 0; i < 2; ++i) + { + auto& p0 = pts[0]; + auto& p1 = pts[i + 1]; + auto& p2 = pts[2 - i]; + auto planeNormal = internal::cross(triangleNormal, p2 - p0); + + component(pcoords, i) = static_cast(internal::dot(wc - p0, planeNormal)) / + static_cast(internal::dot(p1 - p0, planeNormal)); + } + + return ErrorCode::SUCCESS; +} + +} //namespace vtkc + +#endif //vtk_c_Triangle_h diff --git a/vtkc/Vertex.h b/vtkc/Vertex.h new file mode 100644 index 000000000..b06693c51 --- /dev/null +++ b/vtkc/Vertex.h @@ -0,0 +1,131 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.md 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. +//============================================================================ +#ifndef vtk_c_Vertex_h +#define vtk_c_Vertex_h + +#include +#include + +namespace vtkc +{ + +class Vertex : public Cell +{ +public: + constexpr VTKC_EXEC Vertex() : Cell(ShapeId::VERTEX, 1) {} + constexpr VTKC_EXEC explicit Vertex(const Cell& cell) : Cell(cell) {} +}; + +VTKC_EXEC inline vtkc::ErrorCode validate(Vertex tag) noexcept +{ + if (tag.shape() != ShapeId::VERTEX) + { + return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE; + } + if (tag.numberOfPoints() != 1) + { + return ErrorCode::INVALID_NUMBER_OF_POINTS; + } + + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricCenter(Vertex, CoordType&&) noexcept +{ + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricPoint(Vertex, IdComponent pointId, CoordType&&) noexcept +{ + if (pointId == 0) + { + return ErrorCode::SUCCESS; + } + else + { + return ErrorCode::INVALID_POINT_ID; + } +} + +template +VTKC_EXEC inline ComponentType parametricDistance(Vertex, const CoordType&) noexcept +{ + return ComponentType{1}; +} + +template +VTKC_EXEC inline bool cellInside(Vertex, const CoordType&) noexcept +{ + return false; +} + +template +VTKC_EXEC inline vtkc::ErrorCode interpolate( + Vertex, + const Values& values, + const CoordType&, + Result&& result) noexcept +{ + using T = ComponentType; + for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c) + { + component(result, c) = static_cast(values.getValue(0, c)); + } + + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline vtkc::ErrorCode derivative( + Vertex, + const Points&, + const Values& values, + const CoordType&, + Result&& dx, + Result&& dy, + Result&& dz) noexcept +{ + using T = ComponentType; + for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c) + { + component(dx, c) = component(dy, c) = component(dz, c) = T{0}; + } + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricToWorld( + Vertex, + const Points& points, + const PCoordType&, + WCoordType&& wcoords) noexcept +{ + using T = ComponentType; + component(wcoords, 0) = static_cast(points.getValue(0, 0)); + component(wcoords, 1) = static_cast(points.getValue(0, 1)); + component(wcoords, 2) = static_cast(points.getValue(0, 2)); + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline vtkc::ErrorCode worldToParametric( + Vertex, + const Points&, + const WCoordType&, + PCoordType&&) noexcept +{ + return ErrorCode::SUCCESS; +} + +} //namespace vtkc + +#endif //vtk_c_Vertex_h diff --git a/vtkc/Voxel.h b/vtkc/Voxel.h new file mode 100644 index 000000000..b4807b4ff --- /dev/null +++ b/vtkc/Voxel.h @@ -0,0 +1,136 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.md 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. +//============================================================================ +#ifndef vtk_c_Voxel_h +#define vtk_c_Voxel_h + +#include +#include + +#include + +namespace vtkc +{ + +class Voxel : public Hexahedron +{ +public: + constexpr VTKC_EXEC Voxel() : Hexahedron(Cell(ShapeId::VOXEL, 8)) {} + constexpr VTKC_EXEC explicit Voxel(const Cell& cell) : Hexahedron(cell) {} +}; + +namespace internal +{ + +template +VTKC_EXEC inline int getVoxelSpacing(const Points& points, T spacing[3]) +{ + int zeros = 0; + for (int i = 0; i < 3; ++i) + { + spacing[i] = static_cast(points.getValue(6, i) - points.getValue(0, i)); + if (spacing[i] == T(0)) + { + zeros |= 1 << i; + } + } + return zeros; +} + +} // internal + +template +VTKC_EXEC inline vtkc::ErrorCode derivative( + Voxel, + const Points& points, + const Values& values, + const CoordType& pcoords, + Result&& dx, + Result&& dy, + Result&& dz) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + using ProcessingType = internal::ClosestFloatType; + using ResultCompType = ComponentType; + + ProcessingType spacing[3]; + if (internal::getVoxelSpacing(points, spacing) != 0) + { + return ErrorCode::DEGENERATE_CELL_DETECTED; + } + + for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c) + { + ProcessingType dvdp[3]; + internal::parametricDerivative(Hexahedron{}, values, c, pcoords, dvdp); + + component(dx, c) = static_cast(dvdp[0] / spacing[0]); + component(dy, c) = static_cast(dvdp[1] / spacing[1]); + component(dz, c) = static_cast(dvdp[2] / spacing[2]); + } + + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC vtkc::ErrorCode parametricToWorld( + Voxel, + const Points& points, + const PCoordType& pcoords, + WCoordType&& wcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType); + + using T = ComponentType; + + T spacing[3]; + if (internal::getVoxelSpacing(points, spacing) != 0) + { + return ErrorCode::DEGENERATE_CELL_DETECTED; + } + + component(wcoords, 0) = points.getValue(0, 0) + + (spacing[0] * static_cast(component(pcoords, 0))); + component(wcoords, 1) = points.getValue(0, 1) + + (spacing[1] * static_cast(component(pcoords, 1))); + component(wcoords, 2) = points.getValue(0, 2) + + (spacing[2] * static_cast(component(pcoords, 2))); + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC vtkc::ErrorCode worldToParametric( + Voxel, + const Points& points, + const WCoordType& wcoords, + PCoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType); + + using T = ComponentType; + + T spacing[3]; + if (internal::getVoxelSpacing(points, spacing) != 0) + { + return ErrorCode::DEGENERATE_CELL_DETECTED; + } + + component(pcoords, 0) = static_cast(component(wcoords, 0) - points.getValue(0, 0)) / + spacing[0]; + component(pcoords, 1) = static_cast(component(wcoords, 1) - points.getValue(0, 1)) / + spacing[1]; + component(pcoords, 2) = static_cast(component(wcoords, 2) - points.getValue(0, 2)) / + spacing[2]; + return ErrorCode::SUCCESS; +} + +} //namespace vtkc + +#endif //vtk_c_Voxel_h diff --git a/vtkc/Wedge.h b/vtkc/Wedge.h new file mode 100644 index 000000000..03dc3a80c --- /dev/null +++ b/vtkc/Wedge.h @@ -0,0 +1,228 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.md 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. +//============================================================================ +#ifndef vtk_c_Wedge_h +#define vtk_c_Wedge_h + +#include +#include + +#include + +namespace vtkc +{ + +class Wedge : public Cell +{ +public: + constexpr VTKC_EXEC Wedge() : Cell(ShapeId::WEDGE, 6) {} + constexpr VTKC_EXEC explicit Wedge(const Cell& cell) : Cell(cell) {} +}; + +VTKC_EXEC inline vtkc::ErrorCode validate(Wedge tag) noexcept +{ + if (tag.shape() != ShapeId::WEDGE) + { + return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE; + } + if (tag.numberOfPoints() != 6) + { + return ErrorCode::INVALID_NUMBER_OF_POINTS; + } + + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricCenter(Wedge, CoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + using T = ComponentType; + component(pcoords, 0) = T(1)/T(3); + component(pcoords, 1) = T(1)/T(3); + component(pcoords, 2) = 0.5f; + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricPoint( + Wedge, IdComponent pointId, CoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + switch (pointId) + { + case 0: + component(pcoords, 0) = 0.0f; + component(pcoords, 1) = 0.0f; + component(pcoords, 2) = 0.0f; + break; + case 1: + component(pcoords, 0) = 1.0f; + component(pcoords, 1) = 0.0f; + component(pcoords, 2) = 0.0f; + break; + case 2: + component(pcoords, 0) = 0.0f; + component(pcoords, 1) = 1.0f; + component(pcoords, 2) = 0.0f; + break; + case 3: + component(pcoords, 0) = 0.0f; + component(pcoords, 1) = 0.0f; + component(pcoords, 2) = 1.0f; + break; + case 4: + component(pcoords, 0) = 1.0f; + component(pcoords, 1) = 0.0f; + component(pcoords, 2) = 1.0f; + break; + case 5: + component(pcoords, 0) = 0.0f; + component(pcoords, 1) = 1.0f; + component(pcoords, 2) = 1.0f; + break; + default: + return ErrorCode::INVALID_POINT_ID; + } + + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline ComponentType parametricDistance(Wedge, const CoordType& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + return internal::findParametricDistance(pcoords, 3); +} + +template +VTKC_EXEC inline bool cellInside(Wedge, const CoordType& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + using T = ComponentType; + return component(pcoords, 0) >= T{0} && + component(pcoords, 1) >= T{0} && + component(pcoords, 2) >= T{0} && + (component(pcoords, 0) + component(pcoords, 1)) <= T{1} && + component(pcoords, 2) <= T{1}; +} + +template +VTKC_EXEC vtkc::ErrorCode interpolate( + Wedge, + const Values& values, + const CoordType& pcoords, + Result&& result) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + using T = internal::ClosestFloatType; + + auto p0 = static_cast(component(pcoords, 0)); + auto p1 = static_cast(component(pcoords, 1)); + auto p2 = static_cast(component(pcoords, 2)); + auto sm = T{1} - p0 - p1; + + for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c) + { + T t1 = static_cast(values.getValue(0, c)) * sm + + static_cast(values.getValue(1, c)) * p0 + + static_cast(values.getValue(2, c)) * p1; + T t2 = static_cast(values.getValue(3, c)) * sm + + static_cast(values.getValue(4, c)) * p0 + + static_cast(values.getValue(5, c)) * p1; + component(result, c) = static_cast>(internal::lerp(t1, t2, p2)); + } + + return ErrorCode::SUCCESS; +} + +namespace internal +{ + +template +VTKC_EXEC inline void parametricDerivative( + Wedge, const Values& values, IdComponent comp, const CoordType& pcoords, Result&& result) noexcept +{ + using T = internal::ClosestFloatType; + auto p0 = static_cast(component(pcoords, 0)); + auto p1 = static_cast(component(pcoords, 1)); + auto p2 = static_cast(component(pcoords, 2)); + auto rm = T{1} - p2; + auto sm = T{1} - p0 - p1; + + T dr = (static_cast(values.getValue(0, comp)) * -rm) + + (static_cast(values.getValue(1, comp)) * rm) + + (static_cast(values.getValue(3, comp)) * -p2) + + (static_cast(values.getValue(4, comp)) * p2); + + T ds = (static_cast(values.getValue(0, comp)) * -rm) + + (static_cast(values.getValue(2, comp)) * rm) + + (static_cast(values.getValue(3, comp)) * -p2) + + (static_cast(values.getValue(5, comp)) * p2); + + T dt = (static_cast(values.getValue(0, comp)) * -sm) + + (static_cast(values.getValue(1, comp)) * -p0) + + (static_cast(values.getValue(2, comp)) * -p1) + + (static_cast(values.getValue(3, comp)) * sm) + + (static_cast(values.getValue(4, comp)) * p0) + + (static_cast(values.getValue(5, comp)) * p1); + + component(result, 0) = static_cast>(dr); + component(result, 1) = static_cast>(ds); + component(result, 2) = static_cast>(dt); +} + +} // internal + +template +VTKC_EXEC inline vtkc::ErrorCode derivative( + Wedge, + const Points& points, + const Values& values, + const CoordType& pcoords, + Result&& dx, + Result&& dy, + Result&& dz) noexcept +{ + return internal::derivative3D(Wedge{}, + points, + values, + pcoords, + std::forward(dx), + std::forward(dy), + std::forward(dz)); +} + +template +VTKC_EXEC inline vtkc::ErrorCode parametricToWorld( + Wedge, + const Points& points, + const PCoordType& pcoords, + WCoordType&& wcoords) noexcept +{ + return interpolate(Wedge{}, points, pcoords, std::forward(wcoords)); +} + +template +VTKC_EXEC inline vtkc::ErrorCode worldToParametric( + Wedge, + const Points& points, + const WCoordType& wcoords, + PCoordType&& pcoords) noexcept +{ + return internal::worldToParametric3D(Wedge{}, points, wcoords, std::forward(pcoords)); +} + +} // vtkc + +#endif // vtk_c_Wedge_h diff --git a/vtkc/internal/Common.h b/vtkc/internal/Common.h new file mode 100644 index 000000000..85d280939 --- /dev/null +++ b/vtkc/internal/Common.h @@ -0,0 +1,310 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.md 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. +//============================================================================ +#ifndef vtk_c_internal_Common_h +#define vtk_c_internal_Common_h + +#include +#include + +#include +#include + +#define VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType) \ + static_assert(std::is_floating_point>::value, \ + "parametric coordinates should be of floating point type"); + +namespace vtkc +{ +namespace internal +{ + +///========================================================================= +template +VTKC_EXEC inline ComponentType findParametricDistance( + const CoordType& pvals, IdComponent numVals) noexcept +{ + using T = ComponentType; + + T pDistMax = 0.0f; + for (IdComponent i = 0; i < numVals; ++i) + { + ComponentType pDist = 0.0f; + if (component(pvals, i) < T(0)) + { + pDist = -1.0f * component(pvals, i); + } + else if (component(pvals, i) > T(1)) + { + pDist = component(pvals, i) - 1.0f; + } + if (pDist > pDistMax) + { + pDistMax = pDist; + } + } + + return pDistMax; +} + +///========================================================================= +/// Forward declaration +#define FORWARD_DECLAR_PARAMETRIC_DERIVATIVE(tag) \ +template \ +VTKC_EXEC inline void parametricDerivative( \ + vtkc::tag, const Values& values, IdComponent comp, const CoordType&, Result&& result) noexcept + +FORWARD_DECLAR_PARAMETRIC_DERIVATIVE(Triangle); +FORWARD_DECLAR_PARAMETRIC_DERIVATIVE(Quad); +FORWARD_DECLAR_PARAMETRIC_DERIVATIVE(Tetra); +FORWARD_DECLAR_PARAMETRIC_DERIVATIVE(Hexahedron); +FORWARD_DECLAR_PARAMETRIC_DERIVATIVE(Wedge); +FORWARD_DECLAR_PARAMETRIC_DERIVATIVE(Pyramid); + +#undef FORWARD_DECLAR_PARAMETRIC_DERIVATIVE + +///========================================================================= +template +class Space2D +{ +public: + explicit VTKC_EXEC Space2D(const Vector& origin, const Vector& p1, const Vector& p2) noexcept + { + this->Origin = origin; + this->XAxis = p1 - origin; + auto normal = internal::cross(this->XAxis, p2 - origin); + this->YAxis = internal::cross(normal, this->XAxis); + + internal::normalize(this->XAxis); + internal::normalize(this->YAxis); + } + + VTKC_EXEC Vector to2DPoint(Vector pt) const noexcept + { + pt -= this->Origin; + return Vector{ internal::dot(pt, this->XAxis), internal::dot(pt, this->YAxis) }; + } + + VTKC_EXEC Vector to3DVec(const Vector& vec) const noexcept + { + return (this->XAxis * vec[0]) + (this->YAxis * vec[1]); + } + +private: + Vector Origin; + Vector XAxis, YAxis; +}; + +template +VTKC_EXEC inline void jacobian2D( + CellTag tag, const Points& points, const PCoords& pcoords, Matrix& jacobian) noexcept +{ + T pd[2]; + parametricDerivative(tag, points, 0, pcoords, pd); + jacobian(0, 0) = pd[0]; + jacobian(0, 1) = pd[1]; + parametricDerivative(tag, points, 1, pcoords, pd); + jacobian(1, 0) = pd[0]; + jacobian(1, 1) = pd[1]; +} + +template +VTKC_EXEC inline vtkc::ErrorCode derivative2D( + CellTag tag, + const Points& points, + const Values& values, + const CoordType& pcoords, + Result&& dx, + Result&& dy, + Result&& dz) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + using ProcessingType = internal::ClosestFloatType; + using ResultCompType = ComponentType; + + constexpr IdComponent numPoints = CellTag{}.numberOfPoints(); + + Vector pts[numPoints]; + for (int i = 0; i < numPoints; ++i) + { + points.getTuple(i, pts[i]); + } + + // 2-D coordinate system on the cell's plane + Space2D planeSpace(pts[0], pts[1], pts[numPoints - 1]); + Vector pts2d[numPoints]; + for (int i = 0; i < numPoints; ++i) + { + pts2d[i] = planeSpace.to2DPoint(pts[i]); + } + + Matrix jacobian; + jacobian2D(tag, makeFieldAccessorNestedSOA(pts2d, 2), pcoords, jacobian); + Matrix invJacobian; + VTKC_RETURN_ON_ERROR(matrixInverse(jacobian, invJacobian)) + + for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c) + { + Vector dvdp; + parametricDerivative(tag, values, c, pcoords, dvdp); + auto d2D = matrixMultiply(dvdp, invJacobian); + auto d3D = planeSpace.to3DVec(d2D); + + component(dx, c) = static_cast(d3D[0]); + component(dy, c) = static_cast(d3D[1]); + component(dz, c) = static_cast(d3D[2]); + } + + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline vtkc::ErrorCode worldToParametric2D( + CellTag tag, + const Points& points, + const WCoordType& wcoords, + PCoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType); + + using TIn = typename Points::ValueType; + using TOut = ComponentType; + + constexpr IdComponent numPoints = CellTag{}.numberOfPoints(); + + Vector pts[numPoints]; + for (int i = 0; i < numPoints; ++i) + { + points.getTuple(i, pts[i]); + } + + // 2-D coordinate system on the cell's plane + Space2D planeSpace(pts[0], pts[1], pts[numPoints - 1]); + Vector pts2d[numPoints]; + for (int i = 0; i < numPoints; ++i) + { + pts2d[i] = planeSpace.to2DPoint(pts[i]); + } + + auto jacobianEvaluator = [&pts2d](const Vector& pc, Matrix& jacobian) { + jacobian2D(CellTag{}, makeFieldAccessorNestedSOA(pts2d, 2), pc, jacobian); + return ErrorCode::SUCCESS; + }; + + auto functionEvaluator = [&points, &planeSpace](const Vector& pc, Vector& wc) { + Vector wc3(0); + VTKC_RETURN_ON_ERROR(parametricToWorld(CellTag{}, points, pc, wc3)) + wc = planeSpace.to2DPoint(wc3); + + return ErrorCode::SUCCESS; + }; + + Vector wcVec{component(wcoords, 0), component(wcoords, 1), component(wcoords, 2)}; + Vector pcVec; + VTKC_RETURN_ON_ERROR(parametricCenter(tag, pcVec)) + VTKC_RETURN_ON_ERROR(newtonsMethod( + jacobianEvaluator, functionEvaluator, planeSpace.to2DPoint(wcVec), pcVec)) + + component(pcoords, 0) = pcVec[0]; + component(pcoords, 1) = pcVec[1]; + + return ErrorCode::SUCCESS; +} + +///========================================================================= +template +VTKC_EXEC inline void jacobian3D( + CellTag tag, const Points& points, const PCoords& pcoords, Matrix& jacobian) noexcept +{ + T pd[3]; + parametricDerivative(tag, points, 0, pcoords, pd); + jacobian(0, 0) = pd[0]; + jacobian(0, 1) = pd[1]; + jacobian(0, 2) = pd[2]; + parametricDerivative(tag, points, 1, pcoords, pd); + jacobian(1, 0) = pd[0]; + jacobian(1, 1) = pd[1]; + jacobian(1, 2) = pd[2]; + parametricDerivative(tag, points, 2, pcoords, pd); + jacobian(2, 0) = pd[0]; + jacobian(2, 1) = pd[1]; + jacobian(2, 2) = pd[2]; +} + +template +VTKC_EXEC inline vtkc::ErrorCode derivative3D( + CellTag tag, + const Points& points, + const Values& values, + const CoordType& pcoords, + Result&& dx, + Result&& dy, + Result&& dz) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + using ProcessingType = internal::ClosestFloatType; + using ResultCompType = ComponentType; + + Matrix jacobian; + jacobian3D(tag, points, pcoords, jacobian); + Matrix invJacobian; + VTKC_RETURN_ON_ERROR(matrixInverse(jacobian, invJacobian)) + + for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c) + { + Vector dvdp; + parametricDerivative(tag, values, c, pcoords, dvdp); + auto deriv = matrixMultiply(dvdp, invJacobian); + component(dx, c) = static_cast(deriv[0]); + component(dy, c) = static_cast(deriv[1]); + component(dz, c) = static_cast(deriv[2]); + } + + return ErrorCode::SUCCESS; +} + +template +VTKC_EXEC inline vtkc::ErrorCode worldToParametric3D( + CellTag tag, + const Points& points, + const WCoordType& wcoords, + PCoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType); + + using TIn = typename Points::ValueType; + using TOut = ComponentType; + + auto jacobianEvaluator = [tag, &points](const Vector& pc, Matrix& jacobian) { + jacobian3D(tag, points, pc, jacobian); + return ErrorCode::SUCCESS; + }; + + auto functionEvaluator = [tag, &points](const Vector& pc, Vector& wc) { + return parametricToWorld(tag, points, pc, wc); + }; + + internal::Vector wcVec{component(wcoords, 0), component(wcoords, 1), component(wcoords, 2)}; + internal::Vector pcVec; + VTKC_RETURN_ON_ERROR(parametricCenter(tag, pcVec)) + VTKC_RETURN_ON_ERROR(newtonsMethod(jacobianEvaluator, functionEvaluator, wcVec, pcVec)) + + component(pcoords, 0) = pcVec[0]; + component(pcoords, 1) = pcVec[1]; + component(pcoords, 2) = pcVec[2]; + + return ErrorCode::SUCCESS; +} + +} +} // vtkc::internal + +#endif //vtk_c_internal_Common_h diff --git a/vtkc/internal/Config.h b/vtkc/internal/Config.h new file mode 100644 index 000000000..005623cca --- /dev/null +++ b/vtkc/internal/Config.h @@ -0,0 +1,38 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.md 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. +//============================================================================ +#ifndef vtk_c_internal_Config_h +#define vtk_c_internal_Config_h + +#include +#include + +#ifdef __CUDACC__ +# define VTKC_EXEC __device__ __host__ +#else +# define VTKC_EXEC +#endif + +namespace vtkc +{ + +namespace internal +{ +template +using ClosestFloatType = + typename std::enable_if::value, + typename std::conditional::type>::type; +} + +using IdShape = std::int8_t; +using IdComponent = std::int32_t; + +} // namespace vtkc + +#endif // vtk_c_internal_Config_h diff --git a/vtkc/internal/Math.h b/vtkc/internal/Math.h new file mode 100644 index 000000000..99e17a116 --- /dev/null +++ b/vtkc/internal/Math.h @@ -0,0 +1,712 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.md 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. +//============================================================================ +#ifndef vtk_c_internal_Math_h +#define vtk_c_internal_Math_h + +#include + +#ifdef __CUDA_ARCH__ +# define VTKC_MATH_CALL(func, ...) func(__VA_ARGS__) +#else +# define VTKC_MATH_CALL(func, ...) std::func(__VA_ARGS__) +#endif + +namespace vtkc +{ + +namespace internal +{ +///========================================================================= +/// Basic vector class and functions for internal use + +template +class Vector +{ +public: +#if defined(_MSC_VER) && (_MSC_VER == 1900) + // workaround vs2015 bug producing the following error: + // error C2476: ‘constexpr’ constructor does not initialize all members + VTKC_EXEC Vector() noexcept {}; +#else + VTKC_EXEC + constexpr Vector() noexcept {}; +#endif + + VTKC_EXEC + explicit Vector(const T& val) noexcept + { + for (auto& c : Data) + { + c = val; + } + } + + template + VTKC_EXEC + constexpr explicit Vector(const T& c1, const Ts&... cs) noexcept + : Data{c1, cs...} + { + static_assert(sizeof...(Ts) == (Dim - 1), "Invalid number of components passed"); + } + + static constexpr IdComponent getNumberOfComponents() noexcept + { + return Dim; + } + + VTKC_EXEC + T& operator[](IdComponent c) noexcept + { + return this->Data[c]; + } + + VTKC_EXEC + constexpr const T& operator[](IdComponent c) const noexcept + { + return this->Data[c]; + } + + VTKC_EXEC + T* data() noexcept + { + return this->Data; + } + + VTKC_EXEC + constexpr const T* data() const noexcept + { + return this->Data; + } + +private: + T Data[Dim]; +}; + +//--------------------------------------------------------------------------- +template +VTKC_EXEC +Vector& operator+=(Vector& v1, const Vector& v2) noexcept +{ + for (IdComponent i = 0; i < Dim; ++i) + { + v1[i] += v2[i]; + } + return v1; +} + +template +VTKC_EXEC +Vector& operator-=(Vector& v1, const Vector& v2) noexcept +{ + for (IdComponent i = 0; i < Dim; ++i) + { + v1[i] -= v2[i]; + } + return v1; +} + +template +VTKC_EXEC +Vector& operator*=(Vector& v1, const Vector& v2) noexcept +{ + for (IdComponent i = 0; i < Dim; ++i) + { + v1[i] *= v2[i]; + } + return v1; +} + +template +VTKC_EXEC +Vector& operator*=(Vector& v, const T& s) noexcept +{ + for (IdComponent i = 0; i < Dim; ++i) + { + v[i] *= s; + } + return v; +} + +template +VTKC_EXEC +Vector& operator/=(Vector& v1, const Vector& v2) noexcept +{ + for (IdComponent i = 0; i < Dim; ++i) + { + v1[i] /= v2[i]; + } + return v1; +} + +template +VTKC_EXEC +Vector& operator/=(Vector& v, const T& s) noexcept +{ + for (IdComponent i = 0; i < Dim; ++i) + { + v[i] /= s; + } + return v; +} + +template +VTKC_EXEC +Vector operator+(Vector v1, const Vector& v2) noexcept +{ + return v1 += v2; +} + +template +VTKC_EXEC +Vector operator-(Vector v1, const Vector& v2) noexcept +{ + return v1 -= v2; +} + +template +VTKC_EXEC +Vector operator*(Vector v1, const Vector& v2) noexcept +{ + return v1 *= v2; +} + +template +VTKC_EXEC +Vector operator*(Vector v, const T& s) noexcept +{ + return v *= s; +} + +template +VTKC_EXEC +Vector operator*(const T& s, Vector v) noexcept +{ + return v *= s; +} + +template +VTKC_EXEC +Vector operator/(Vector v1, const Vector& v2) noexcept +{ + return v1 /= v2; +} + +template +VTKC_EXEC +Vector operator/(Vector v, const T& s) noexcept +{ + return v /= s; +} + +//--------------------------------------------------------------------------- +template +VTKC_EXEC +T dot(const Vector& v1, const Vector& v2) noexcept +{ + T result{}; + for (IdComponent i = 0; i < Dim; ++i) + { + result += v1[i] * v2[i]; + } + return result; +} + +template +VTKC_EXEC +Vector cross(const Vector& v1, const Vector& v2) noexcept +{ + return Vector(v1[1] * v2[2] - v1[2] * v2[1], + v1[2] * v2[0] - v1[0] * v2[2], + v1[0] * v2[1] - v1[1] * v2[0]); +} + +template +VTKC_EXEC +T magnitude(const Vector& v) +{ + return VTKC_MATH_CALL(sqrt, (dot(v, v))); +} + +template +VTKC_EXEC +void normalize(Vector& v) +{ + v /= magnitude(v); +} + +template +VTKC_EXEC +Vector normal(Vector v) +{ + normalize(v); + return v; +} + +///========================================================================= +/// Basic matrix class and functions for internal use + +template +class Matrix +{ +public: +#if defined(_MSC_VER) && (_MSC_VER == 1900) + // workaround vs2015 bug producing the following error: + // error C2476: ‘constexpr’ constructor does not initialize all members + VTKC_EXEC Matrix() noexcept {}; +#else + VTKC_EXEC + constexpr Matrix() noexcept {}; +#endif + + VTKC_EXEC + T& operator()(IdComponent r, IdComponent c) noexcept + { + return this->Columns[c][r]; + } + + VTKC_EXEC + constexpr const T& operator()(IdComponent r, IdComponent c) const noexcept + { + return this->Columns[c][r]; + } + + VTKC_EXEC + static constexpr IdComponent getNumberOfRows() noexcept + { + return NumberOfRows; + } + + VTKC_EXEC + static constexpr IdComponent getNumberOfColumns() noexcept + { + return NumberOfColumns; + } + + VTKC_EXEC + constexpr const Vector& getColumn(IdComponent c) const + { + return this->Columns[c]; + } + + VTKC_EXEC + Vector getRow(IdComponent r) const + { + Vector row; + for (IdComponent i = 0; i < NumberOfColumns; ++i) + { + row[i] = this->Columns[i][r]; + } + return row; + } + + VTKC_EXEC + Matrix& operator+=(const Matrix& m) noexcept + { + for (IdComponent i = 0; i < NumberOfColumns; ++i) + { + this->Columns[i] += m.Columns[i]; + } + return *this; + } + + VTKC_EXEC + Matrix& operator-=(const Matrix& m) noexcept + { + for (IdComponent i = 0; i < NumberOfColumns; ++i) + { + this->Columns[i] -= m.Columns[i]; + } + return *this; + } + + VTKC_EXEC + Matrix& operator*=(const T& s) noexcept + { + for (auto& c : this->Columns) + { + c *= s; + } + return *this; + } + + VTKC_EXEC + Matrix& operator/=(const T& s) noexcept + { + for (auto& c : this->Columns) + { + c /= s; + } + return *this; + } + +private: + Vector Columns[NumberOfColumns]; +}; + +//--------------------------------------------------------------------------- +template +VTKC_EXEC +Matrix operator+(Matrix m1, + const Matrix& m2) noexcept +{ + return m1 += m2; +} + +template +VTKC_EXEC +Matrix operator-(Matrix m1, + const Matrix& m2) noexcept +{ + return m1 -= m2; +} + +template +VTKC_EXEC +Matrix operator*(Matrix m, const T& s) noexcept +{ + return m *= s; +} + +template +VTKC_EXEC +Matrix operator/(Matrix m, const T& s) noexcept +{ + return m /= s; +} + +//--------------------------------------------------------------------------- +template +VTKC_EXEC +Matrix matrixMultiply( + const Matrix& leftFactor, + const Matrix& rightFactor) noexcept +{ + Matrix result; + for (IdComponent rowIndex = 0; rowIndex < NumRow; ++rowIndex) + { + for (IdComponent colIndex = 0; colIndex < NumCol; ++colIndex) + { + T sum = T(leftFactor(rowIndex, 0) * rightFactor(0, colIndex)); + for (IdComponent internalIndex = 1; internalIndex < NumInternal; internalIndex++) + { + sum = T(sum + (leftFactor(rowIndex, internalIndex) * rightFactor(internalIndex, colIndex))); + } + result(rowIndex, colIndex) = sum; + } + } + return result; +} + +template +VTKC_EXEC +Vector matrixMultiply( + const Matrix& leftFactor, + const Vector& rightFactor) noexcept +{ + Vector result; + for (IdComponent rowIndex = 0; rowIndex < NumRow; ++rowIndex) + { + result[rowIndex] = dot(leftFactor.getRow(rowIndex), rightFactor); + } + return result; +} + +template +VTKC_EXEC +Vector matrixMultiply( + const Vector& leftFactor, + const Matrix& rightFactor) noexcept +{ + Vector result; + for (IdComponent colIndex = 0; colIndex < NumCol; ++colIndex) + { + result[colIndex] = dot(leftFactor, rightFactor.getColumn(colIndex)); + } + return result; +} + +//--------------------------------------------------------------------------- +// The following code is copied from VTK-m + +namespace detail +{ + +template +VTKC_EXEC +vtkc::ErrorCode matrixLUPFactorFindPivot(Matrix& m, + Vector& permutation, + IdComponent topCornerIndex, + T& inversionParity) noexcept +{ + constexpr T epsilon = std::is_floating_point::value ? + (std::is_same::value ? T(1e-5f) : T(1e-9)) : + T{}; + + IdComponent maxRowIndex = topCornerIndex; + T maxValue = VTKC_MATH_CALL(abs, (m(maxRowIndex, topCornerIndex))); + for (IdComponent rowIndex = topCornerIndex + 1; rowIndex < Size; rowIndex++) + { + T compareValue = VTKC_MATH_CALL(abs, (m(rowIndex, topCornerIndex))); + if (maxValue < compareValue) + { + maxValue = compareValue; + maxRowIndex = rowIndex; + } + } + + if (maxValue < epsilon) + { + return vtkc::ErrorCode::MATRIX_LUP_FACTORIZATION_FAILED; + } + + if (maxRowIndex != topCornerIndex) + { + // Swap rows in matrix. + for (IdComponent i = 0; i < Size; ++i) + { + auto temp = m(topCornerIndex, i); + m(topCornerIndex, i) = m(maxRowIndex, i); + m(maxRowIndex, i) = temp; + } + + // Record change in permutation matrix. + IdComponent maxOriginalRowIndex = permutation[maxRowIndex]; + permutation[maxRowIndex] = permutation[topCornerIndex]; + permutation[topCornerIndex] = maxOriginalRowIndex; + + // Keep track of inversion parity. + inversionParity = -inversionParity; + } + + return vtkc::ErrorCode::SUCCESS; +} + +// Used with MatrixLUPFactor +template +VTKC_EXEC +void matrixLUPFactorFindUpperTriangleElements( + Matrix& m, IdComponent topCornerIndex) noexcept +{ + // Compute values for upper triangle on row topCornerIndex + for (IdComponent colIndex = topCornerIndex + 1; colIndex < Size; colIndex++) + { + m(topCornerIndex, colIndex) /= m(topCornerIndex, topCornerIndex); + } + + // Update the rest of the matrix for calculations on subsequent rows + for (IdComponent rowIndex = topCornerIndex + 1; rowIndex < Size; rowIndex++) + { + for (IdComponent colIndex = topCornerIndex + 1; colIndex < Size; colIndex++) + { + m(rowIndex, colIndex) -= m(rowIndex, topCornerIndex) * m(topCornerIndex, colIndex); + } + } +} + +/// Performs an LUP-factorization on the given matrix using Crout's method. The +/// LU-factorization takes a matrix A and decomposes it into a lower triangular +/// matrix L and upper triangular matrix U such that A = LU. The +/// LUP-factorization also allows permutation of A, which makes the +/// decomposition always possible so long as A is not singular. In addition to +/// matrices L and U, LUP also finds permutation matrix P containing all zeros +/// except one 1 per row and column such that PA = LU. +/// +/// The result is done in place such that the lower triangular matrix, L, is +/// stored in the lower-left triangle of A including the diagonal. The upper +/// triangular matrix, U, is stored in the upper-right triangle of L not +/// including the diagonal. The diagonal of U in Crout's method is all 1's (and +/// therefore not explicitly stored). +/// +/// The permutation matrix P is represented by the permutation vector. If +/// permutation[i] = j then row j in the original matrix A has been moved to +/// row i in the resulting matrices. The permutation matrix P can be +/// represented by a matrix with p_i,j = 1 if permutation[i] = j and 0 +/// otherwise. If using LUP-factorization to compute a determinant, you also +/// need to know the parity (whether there is an odd or even amount) of +/// inversions. An inversion is an instance of a smaller number appearing after +/// a larger number in permutation. Although you can compute the inversion +/// parity after the fact, this function keeps track of it with much less +/// compute resources. The parameter inversionParity is set to 1.0 for even +/// parity and -1.0 for odd parity. +/// +/// Not all matrices (specifically singular matrices) have an +/// LUP-factorization. If the LUP-factorization succeeds, valid is set to true. +/// Otherwise, valid is set to false and the result is indeterminant. +/// +template +VTKC_EXEC +vtkc::ErrorCode matrixLUPFactor(Matrix& m, + Vector& permutation, + T& inversionParity) noexcept +{ + // Initialize permutation. + for (IdComponent index = 0; index < Size; index++) + { + permutation[index] = index; + } + inversionParity = T(1); + + for (IdComponent rowIndex = 0; rowIndex < Size; rowIndex++) + { + VTKC_RETURN_ON_ERROR(matrixLUPFactorFindPivot(m, permutation, rowIndex, inversionParity)) + matrixLUPFactorFindUpperTriangleElements(m, rowIndex); + } + + return vtkc::ErrorCode::SUCCESS; +} + +/// Use a previous factorization done with MatrixLUPFactor to solve the +/// system Ax = b. Instead of A, this method takes in the LU and P +/// matrices calculated by MatrixLUPFactor from A. The x matrix is returned. +/// +template +VTKC_EXEC +Vector matrixLUPSolve( + const Matrix& LU, + const Vector& permutation, + const Vector& b) noexcept +{ + // The LUP-factorization gives us PA = LU or equivalently A = inv(P)LU. + // Substituting into Ax = b gives us inv(P)LUx = b or LUx = Pb. + // Now consider the intermediate vector y = Ux. + // Substituting in the previous two equations yields Ly = Pb. + // Solving Ly = Pb is easy because L is triangular and P is just a + // permutation. + Vector y; + for (IdComponent rowIndex = 0; rowIndex < Size; rowIndex++) + { + y[rowIndex] = b[permutation[rowIndex]]; + // Recall that L is stored in the lower triangle of LU including diagonal. + for (IdComponent colIndex = 0; colIndex < rowIndex; colIndex++) + { + y[rowIndex] -= LU(rowIndex, colIndex) * y[colIndex]; + } + y[rowIndex] /= LU(rowIndex, rowIndex); + } + + // Now that we have y, we can easily solve Ux = y for x. + Vector x; + for (IdComponent rowIndex = Size - 1; rowIndex >= 0; rowIndex--) + { + // Recall that U is stored in the upper triangle of LU with the diagonal + // implicitly all 1's. + x[rowIndex] = y[rowIndex]; + for (IdComponent colIndex = rowIndex + 1; colIndex < Size; colIndex++) + { + x[rowIndex] -= LU(rowIndex, colIndex) * x[colIndex]; + } + } + + return x; +} + +} // namespace detail + +template +VTKC_EXEC +vtkc::ErrorCode solveLinearSystem(const Matrix& A, + const Vector& b, + Vector& x) noexcept +{ + // First, we will make an LUP-factorization to help us. + Matrix LU = A; + Vector permutation; + T inversionParity; // Unused. + VTKC_RETURN_ON_ERROR(detail::matrixLUPFactor(LU, permutation, inversionParity)) + + // Next, use the decomposition to solve the system. + x = detail::matrixLUPSolve(LU, permutation, b); + return vtkc::ErrorCode::SUCCESS; +} + +/// Find and return the inverse of the given matrix. If the matrix is singular, +/// the inverse will not be correct and valid will be set to false. +/// +template +VTKC_EXEC +vtkc::ErrorCode matrixInverse(const Matrix& A, Matrix& invA) noexcept +{ + // First, we will make an LUP-factorization to help us. + Matrix LU = A; + Vector permutation; + T inversionParity; // Unused + VTKC_RETURN_ON_ERROR(detail::matrixLUPFactor(LU, permutation, inversionParity)) + + // We will use the decomposition to solve AX = I for X where X is + // clearly the inverse of A. Our solve method only works for vectors, + // so we solve for one column of invA at a time. + Vector ICol(T(0)); + for (IdComponent colIndex = 0; colIndex < Size; colIndex++) + { + ICol[colIndex] = 1; + Vector invACol = detail::matrixLUPSolve(LU, permutation, ICol); + ICol[colIndex] = 0; + for (IdComponent i = 0; i < Size; ++i) + { + invA(i, colIndex) = invACol[i]; + } + } + return vtkc::ErrorCode::SUCCESS; +} + +///========================================================================= +template +VTKC_EXEC +inline vtkc::ErrorCode newtonsMethod( + const JacobianFunctor& jacobianEvaluator, + const FunctionFunctor& functionEvaluator, + const Vector& rhs, + Vector& result, + T convergeDifference = 1e-3f, + int maxIterations = 10) noexcept +{ + auto x = result; // pass initial guess in result + bool converged = false; + + for (int i = 0; !converged && i < maxIterations; ++i) + { + Matrix jacobian; + Vector fx(0); + + VTKC_RETURN_ON_ERROR(jacobianEvaluator(x, jacobian)) + VTKC_RETURN_ON_ERROR(functionEvaluator(x, fx)) + + Vector deltax; + VTKC_RETURN_ON_ERROR(solveLinearSystem(jacobian, fx - rhs, deltax)) + x -= deltax; + + converged = true; + for (int c = 0; c < Size; ++c) + { + converged &= (VTKC_MATH_CALL(abs, (deltax[c])) < convergeDifference); + } + } + + result = x; + return converged ? ErrorCode::SUCCESS : ErrorCode::SOLUTION_DID_NOT_CONVERGE; +} + +///========================================================================= +template +VTKC_EXEC inline T lerp(T v0, T v1, T t) +{ + static_assert(std::is_floating_point::value, + "lerp requires floating point parameters"); + + return VTKC_MATH_CALL(fma, (t), (v1), (VTKC_MATH_CALL(fma, (-t), (v0), (v0)))); +} + +} +} // vtkc::internal + +#endif // vtk_c_internal_Math_h diff --git a/vtkc/vtkc.h b/vtkc/vtkc.h new file mode 100644 index 000000000..61835cf2f --- /dev/null +++ b/vtkc/vtkc.h @@ -0,0 +1,251 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.md 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. +//============================================================================ +#ifndef vtk_c_vtkc_h +#define vtk_c_vtkc_h + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace vtkc +{ + +/// \brief Perform basic checks to validate cell's state. +/// \param[in] tag The cell tag to validate. +/// \return vtkc::ErrorCode::SUCCESS if valid. +/// +VTKC_EXEC inline vtkc::ErrorCode validate(Cell tag) noexcept +{ + ErrorCode status = ErrorCode::SUCCESS; + switch (tag.shape()) + { + vtkcGenericCellShapeMacro(status = validate(CellTag{tag})); + default: + status = ErrorCode::INVALID_SHAPE_ID; + break; + } + return status; +} + +/// \brief Return center of a cell in parametric coordinates. +/// \remark Note that the parametric center is not always located at (0.5,0.5,0.5). +/// \param[in] tag The cell tag. +/// \param[out] pcoords The center of the cell in parametric coordinates. +/// \return A vtkc::ErrorCode value indicating the status of the operation. +/// +template +VTKC_EXEC inline vtkc::ErrorCode parametricCenter(Cell tag, CoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + ErrorCode status = ErrorCode::SUCCESS; + switch (tag.shape()) + { + vtkcGenericCellShapeMacro(status = parametricCenter(CellTag{tag}, std::forward(pcoords))); + default: + status = ErrorCode::INVALID_SHAPE_ID; + break; + } + return status; +} + +/// \brief Return the parametric coordinates of a cell's point. +/// \param[in] tag The cell tag. +/// \param[in] pointId The point number. +/// \param[out] pcoords The parametric coordinates of a cell's point. +/// \return A vtkc::ErrorCode value indicating the status of the operation. +/// +template +VTKC_EXEC inline vtkc::ErrorCode parametricPoint( + Cell tag, IdComponent pointId, CoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + ErrorCode status = ErrorCode::SUCCESS; + switch (tag.shape()) + { + vtkcGenericCellShapeMacro(status = parametricPoint(CellTag{tag}, pointId, std::forward(pcoords))); + default: + status = ErrorCode::INVALID_SHAPE_ID; + break; + } + return status; +} + +/// \brief Return the parametric distance of a parametric coordinate to a cell. +/// \param[in] tag The cell tag. +/// \param[in] pcoords The parametric coordinates of the point. +/// \return The parametric distance of the point to the cell. +/// If point is inside the cell, 0 is returned. +/// \pre tag should be a valid cell, otherwise the result is undefined. +/// +template +VTKC_EXEC inline ComponentType parametricDistance(Cell tag, const CoordType& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + ComponentType dist{0}; + switch (tag.shape()) + { + vtkcGenericCellShapeMacro(dist = parametricDistance(CellTag{tag}, pcoords)); + default: + break; + } + return dist; +} + +/// \brief Check if the given parametric point lies inside a cell. +/// \param[in] tag The cell tag. +/// \param[in] pcoords The parametric coordinates of the point. +/// \return true if inside, false otherwise. +/// \pre tag should be a valid cell, otherwise the result is undefined. +/// +template +VTKC_EXEC inline bool cellInside(Cell tag, const CoordType& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + bool inside = false; + switch (tag.shape()) + { + vtkcGenericCellShapeMacro(inside = cellInside(CellTag{tag}, pcoords)); + default: + break; + } + return inside; +} + +/// \brief Interpolate \c values at the paramteric coordinates \c pcoords +/// \param[in] tag The cell tag. +/// \param[in] values A \c FieldAccessor for values to interpolate +/// \param[in] pcoords The parametric coordinates. +/// \param[out] result The interpolation result +/// \return A vtkc::ErrorCode value indicating the status of the operation. +/// +template +VTKC_EXEC inline vtkc::ErrorCode interpolate( + Cell tag, + const Values& values, + const CoordType& pcoords, + Result&& result) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + ErrorCode status = ErrorCode::SUCCESS; + switch (tag.shape()) + { + vtkcGenericCellShapeMacro(status = interpolate(CellTag{tag}, values, pcoords, std::forward(result))); + default: + status = ErrorCode::INVALID_SHAPE_ID; + } + return status; +} + +/// \brief Compute derivative of \c values at the paramteric coordinates \c pcoords +/// \param[in] tag The cell tag. +/// \param[in] points A \c FieldAccessor for points of the cell +/// \param[in] values A \c FieldAccessor for the values to compute derivative of. +/// \param[in] pcoords The parametric coordinates. +/// \param[out] dx The derivative along X +/// \param[out] dy The derivative along Y +/// \param[out] dz The derivative along Z +/// \return A vtkc::ErrorCode value indicating the status of the operation. +/// +template +VTKC_EXEC inline vtkc::ErrorCode derivative( + Cell tag, + const Points& points, + const Values& values, + const CoordType& pcoords, + Result&& dx, + Result&& dy, + Result&& dz) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType); + + ErrorCode status = ErrorCode::SUCCESS; + switch (tag.shape()) + { + vtkcGenericCellShapeMacro(status = derivative(CellTag{tag}, + points, + values, + pcoords, + std::forward(dx), + std::forward(dy), + std::forward(dz))); + default: + status = ErrorCode::INVALID_SHAPE_ID; + } + return status; +} + +/// \brief Compute world coordinates from parametric coordinates +/// \param[in] tag The cell tag. +/// \param[in] points A \c FieldAccessor for points of the cell +/// \param[in] pcoords The parametric coordinates. +/// \param[out] wcoords The world coordinates. +/// +template +VTKC_EXEC inline vtkc::ErrorCode parametricToWorld( + Cell tag, + const Points& points, + const PCoordType& pcoords, + WCoordType&& wcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType); + + ErrorCode status = ErrorCode::SUCCESS; + switch (tag.shape()) + { + vtkcGenericCellShapeMacro(status = parametricToWorld(CellTag{tag}, points, pcoords, std::forward(wcoords))); + default: + status = ErrorCode::INVALID_SHAPE_ID; + } + return status; +} + +/// \brief Compute parametric coordinates from world coordinates +/// \param[in] tag The cell tag. +/// \param[in] points A \c FieldAccessor for points of the cell +/// \param[in] wcoords The world coordinates. +/// \param[out] pcoords The parametric coordinates. +/// +template +VTKC_EXEC inline vtkc::ErrorCode worldToParametric( + Cell tag, + const Points& points, + const WCoordType& wcoords, + PCoordType&& pcoords) noexcept +{ + VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType); + + ErrorCode status = ErrorCode::SUCCESS; + switch (tag.shape()) + { + vtkcGenericCellShapeMacro(status = worldToParametric(CellTag{tag}, points, wcoords, std::forward(pcoords))); + default: + status = ErrorCode::INVALID_SHAPE_ID; + } + return status; +} + +} //namespace vtkc + +#endif //vtk_c_vtkc_h