mirror of
https://gitlab.kitware.com/vtk/vtk-m
synced 2024-09-19 18:45:43 +00:00
vtkc 2019-09-24 (cdc72582)
Code extracted from: https://gitlab.kitware.com/sujin.philip/vtk-c.git at commit cdc72582dae24cc66b3525961e32edbc2a395341 (master).
This commit is contained in:
commit
e909edd61d
43
LICENSE.md
Normal file
43
LICENSE.md
Normal file
@ -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.
|
65
README.md
Normal file
65
README.md
Normal file
@ -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 <vtkc/vtkc.h>
|
||||||
|
|
||||||
|
std::array<float, 3> pcoords;
|
||||||
|
auto status = vtkc::parametricCenter(vtkc::Hexahedron{}, pcoords);
|
||||||
|
|
||||||
|
std::vector<std::array<float, 3>> points = { ... };
|
||||||
|
std::vector<std::array<double, 4>> field = { ... };
|
||||||
|
std::array<double, 4> 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
|
69
vtkc/ErrorCode.h
Normal file
69
vtkc/ErrorCode.h
Normal file
@ -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 <vtkc/internal/Config.h>
|
||||||
|
|
||||||
|
#include <cstdint>
|
||||||
|
|
||||||
|
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
|
367
vtkc/FieldAccessor.h
Normal file
367
vtkc/FieldAccessor.h
Normal file
@ -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 <vtkc/internal/Config.h>
|
||||||
|
|
||||||
|
#include <utility>
|
||||||
|
|
||||||
|
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 <typename T> inline decltype(std::declval<T>()[0]) VecCheck(T);
|
||||||
|
inline NotVecTag VecCheck(...);
|
||||||
|
|
||||||
|
template <typename T>
|
||||||
|
using IsVecType =
|
||||||
|
std::integral_constant<bool, !std::is_same<decltype(VecCheck(std::declval<T>())),
|
||||||
|
internal::NotVecTag>::value>;
|
||||||
|
|
||||||
|
//----------------------------------------------------------------------------
|
||||||
|
/// Does T support operator()(integer_type)?
|
||||||
|
|
||||||
|
struct NotFuncTag {};
|
||||||
|
|
||||||
|
template <typename T> inline decltype(std::declval<T>()(0)) FuncCheck(T);
|
||||||
|
inline NotFuncTag FuncCheck(...);
|
||||||
|
|
||||||
|
template <typename T>
|
||||||
|
using IsFuncType =
|
||||||
|
std::integral_constant<bool, !std::is_same<decltype(FuncCheck(std::declval<T>())),
|
||||||
|
internal::NotFuncTag>::value>;
|
||||||
|
|
||||||
|
//----------------------------------------------------------------------------
|
||||||
|
template <typename T> typename T::size_type IndexTypeImpl(T);
|
||||||
|
IdComponent IndexTypeImpl(...);
|
||||||
|
|
||||||
|
template <typename VecT>
|
||||||
|
using IndexType = decltype(IndexTypeImpl(std::declval<VecT>()));
|
||||||
|
|
||||||
|
//----------------------------------------------------------------------------
|
||||||
|
struct VecTypeIndexer {};
|
||||||
|
struct VecTypeFunctor {};
|
||||||
|
struct VecTypeScalar {};
|
||||||
|
|
||||||
|
template <typename T>
|
||||||
|
using VecTypeTag =
|
||||||
|
typename std::conditional<IsVecType<T>::value,
|
||||||
|
VecTypeIndexer,
|
||||||
|
typename std::conditional<IsFuncType<T>::value,
|
||||||
|
VecTypeFunctor,
|
||||||
|
VecTypeScalar>::type>::type;
|
||||||
|
|
||||||
|
//----------------------------------------------------------------------------
|
||||||
|
template <typename T>
|
||||||
|
VTKC_EXEC auto componentImpl(VecTypeIndexer, T&& vec, int idx) -> decltype(vec[0])
|
||||||
|
{
|
||||||
|
return vec[static_cast<IndexType<T>>(idx)];
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T>
|
||||||
|
VTKC_EXEC auto componentImpl(VecTypeFunctor, T&& vec, int idx) -> decltype(vec(0))
|
||||||
|
{
|
||||||
|
return vec(idx);
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T>
|
||||||
|
VTKC_EXEC T& componentImpl(VecTypeScalar, T&& vec, int)
|
||||||
|
{
|
||||||
|
return vec;
|
||||||
|
}
|
||||||
|
|
||||||
|
} // namespace internal
|
||||||
|
|
||||||
|
//----------------------------------------------------------------------------
|
||||||
|
template <typename T>
|
||||||
|
VTKC_EXEC auto component(T&& vec, int idx)
|
||||||
|
-> decltype(internal::componentImpl<T>(internal::VecTypeTag<T>{}, std::forward<T>(vec), 0))
|
||||||
|
{
|
||||||
|
return internal::componentImpl<T>(internal::VecTypeTag<T>{}, std::forward<T>(vec), idx);
|
||||||
|
}
|
||||||
|
|
||||||
|
//----------------------------------------------------------------------------
|
||||||
|
template <typename T>
|
||||||
|
using ComponentType = typename std::decay<decltype(component(std::declval<T>(), 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 <typename FieldType>
|
||||||
|
/// 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 <typename VecType>
|
||||||
|
/// 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 <typename VecType>
|
||||||
|
/// 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 <typename FieldType>
|
||||||
|
class FieldAccessorNestedSOA
|
||||||
|
{
|
||||||
|
public:
|
||||||
|
using ValueType = typename std::decay<decltype(component(std::declval<FieldType>()[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<IdxType>(tuple)], comp) = value;
|
||||||
|
}
|
||||||
|
|
||||||
|
VTKC_EXEC ValueType getValue(int tuple, int comp) const
|
||||||
|
{
|
||||||
|
return component((*this->Field)[static_cast<IdxType>(tuple)], comp);
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename VecType>
|
||||||
|
VTKC_EXEC void setTuple(int tuple, const VecType& value) const
|
||||||
|
{
|
||||||
|
for (int i = 0; i < this->NumberOfComponents; ++i)
|
||||||
|
{
|
||||||
|
component((*this->Field)[static_cast<IdxType>(tuple)], i) =
|
||||||
|
static_cast<ValueType>(component(value, i));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename VecType>
|
||||||
|
VTKC_EXEC void getTuple(int tuple, VecType& value) const
|
||||||
|
{
|
||||||
|
for (int i = 0; i < this->NumberOfComponents; ++i)
|
||||||
|
{
|
||||||
|
component(value, i) = static_cast<ComponentType<VecType>>(
|
||||||
|
component((*this->Field)[static_cast<IdxType>(tuple)], i));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
private:
|
||||||
|
using IdxType = internal::IndexType<FieldType>;
|
||||||
|
|
||||||
|
FieldType* Field;
|
||||||
|
int NumberOfComponents;
|
||||||
|
};
|
||||||
|
|
||||||
|
template <typename FieldType>
|
||||||
|
class FieldAccessorNestedSOAConst
|
||||||
|
{
|
||||||
|
public:
|
||||||
|
using ValueType = typename std::decay<decltype(component(std::declval<FieldType>()[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<IdxType>(tuple)], comp);
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename VecType>
|
||||||
|
VTKC_EXEC void getTuple(int tuple, VecType& value) const
|
||||||
|
{
|
||||||
|
for (int i = 0; i < this->NumberOfComponents; ++i)
|
||||||
|
{
|
||||||
|
component(value, i) = static_cast<ComponentType<VecType>>(
|
||||||
|
component((*this->Field)[static_cast<IdxType>(tuple)], i));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
private:
|
||||||
|
using IdxType = internal::IndexType<FieldType>;
|
||||||
|
|
||||||
|
const FieldType* Field;
|
||||||
|
int NumberOfComponents;
|
||||||
|
};
|
||||||
|
|
||||||
|
template <typename FieldType>
|
||||||
|
VTKC_EXEC
|
||||||
|
FieldAccessorNestedSOA<FieldType> makeFieldAccessorNestedSOA(FieldType& field,
|
||||||
|
int numberOfComponents = 1)
|
||||||
|
{
|
||||||
|
return FieldAccessorNestedSOA<FieldType>(field, numberOfComponents);
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename FieldType>
|
||||||
|
VTKC_EXEC
|
||||||
|
FieldAccessorNestedSOAConst<FieldType> makeFieldAccessorNestedSOAConst(const FieldType& field,
|
||||||
|
int numberOfComponents = 1)
|
||||||
|
{
|
||||||
|
return FieldAccessorNestedSOAConst<FieldType>(field, numberOfComponents);
|
||||||
|
}
|
||||||
|
|
||||||
|
///----------------------------------------------------------------------------
|
||||||
|
template <typename FieldType>
|
||||||
|
class FieldAccessorFlatSOA
|
||||||
|
{
|
||||||
|
public:
|
||||||
|
using ValueType = typename std::decay<decltype(std::declval<FieldType>()[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<IdxType>(tuple * this->NumberOfComponents + comp);
|
||||||
|
(*this->Field)[FlatIdx] = value;
|
||||||
|
}
|
||||||
|
|
||||||
|
VTKC_EXEC ValueType getValue(int tuple, int comp) const
|
||||||
|
{
|
||||||
|
auto FlatIdx = static_cast<IdxType>(tuple * this->NumberOfComponents + comp);
|
||||||
|
return (*this->Field)[FlatIdx];
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename VecType>
|
||||||
|
VTKC_EXEC void setTuple(int tuple, const VecType& value) const
|
||||||
|
{
|
||||||
|
auto start = static_cast<IdxType>(tuple * this->NumberOfComponents);
|
||||||
|
for (int i = 0; i < this->NumberOfComponents; ++i)
|
||||||
|
{
|
||||||
|
(*this->Field)[start++] = static_cast<ValueType>(component(value, i));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename VecType>
|
||||||
|
VTKC_EXEC void getTuple(int tuple, VecType& value) const
|
||||||
|
{
|
||||||
|
auto start = static_cast<IdxType>(tuple * this->NumberOfComponents);
|
||||||
|
for (int i = 0; i < this->NumberOfComponents; ++i)
|
||||||
|
{
|
||||||
|
component(value, i) = static_cast<ComponentType<VecType>>((*this->Field)[start++]);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
private:
|
||||||
|
using IdxType = internal::IndexType<FieldType>;
|
||||||
|
|
||||||
|
FieldType* Field;
|
||||||
|
int NumberOfComponents;
|
||||||
|
};
|
||||||
|
|
||||||
|
template <typename FieldType>
|
||||||
|
class FieldAccessorFlatSOAConst
|
||||||
|
{
|
||||||
|
public:
|
||||||
|
using ValueType = typename std::decay<decltype(std::declval<FieldType>()[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<IdxType>(tuple * this->NumberOfComponents + comp);
|
||||||
|
return (*this->Field)[FlatIdx];
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename VecType>
|
||||||
|
VTKC_EXEC void getTuple(int tuple, VecType& value) const
|
||||||
|
{
|
||||||
|
auto start = static_cast<IdxType>(tuple * this->NumberOfComponents);
|
||||||
|
for (int i = 0; i < this->NumberOfComponents; ++i)
|
||||||
|
{
|
||||||
|
component(value, i) = static_cast<ComponentType<VecType>>((*this->Field)[start++]);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
private:
|
||||||
|
using IdxType = internal::IndexType<FieldType>;
|
||||||
|
|
||||||
|
const FieldType* Field;
|
||||||
|
int NumberOfComponents;
|
||||||
|
};
|
||||||
|
|
||||||
|
template <typename FieldType>
|
||||||
|
VTKC_EXEC
|
||||||
|
FieldAccessorFlatSOA<FieldType> makeFieldAccessorFlatSOA(FieldType& field,
|
||||||
|
int numberOfComponents)
|
||||||
|
{
|
||||||
|
return FieldAccessorFlatSOA<FieldType>(field, numberOfComponents);
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename FieldType>
|
||||||
|
VTKC_EXEC
|
||||||
|
FieldAccessorFlatSOAConst<FieldType> makeFieldAccessorFlatSOAConst(const FieldType& field,
|
||||||
|
int numberOfComponents)
|
||||||
|
{
|
||||||
|
return FieldAccessorFlatSOAConst<FieldType>(field, numberOfComponents);
|
||||||
|
}
|
||||||
|
|
||||||
|
} // namespace vtkc
|
||||||
|
|
||||||
|
#endif // vtk_c_FieldAccessor_h
|
254
vtkc/Hexahedron.h
Normal file
254
vtkc/Hexahedron.h
Normal file
@ -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 <vtkc/ErrorCode.h>
|
||||||
|
#include <vtkc/Shapes.h>
|
||||||
|
|
||||||
|
#include <vtkc/internal/Common.h>
|
||||||
|
|
||||||
|
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<typename CoordType>
|
||||||
|
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<typename CoordType>
|
||||||
|
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<typename CoordType>
|
||||||
|
VTKC_EXEC inline ComponentType<CoordType> parametricDistance(Hexahedron, const CoordType& pcoords) noexcept
|
||||||
|
{
|
||||||
|
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
|
||||||
|
return internal::findParametricDistance(pcoords, 3);
|
||||||
|
}
|
||||||
|
|
||||||
|
template<typename CoordType>
|
||||||
|
VTKC_EXEC inline bool cellInside(Hexahedron, const CoordType& pcoords) noexcept
|
||||||
|
{
|
||||||
|
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
|
||||||
|
|
||||||
|
using T = ComponentType<CoordType>;
|
||||||
|
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 <typename Values, typename CoordType, typename Result>
|
||||||
|
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<typename Values::ValueType>;
|
||||||
|
|
||||||
|
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
|
||||||
|
{
|
||||||
|
auto vbf = internal::lerp(static_cast<T>(values.getValue(0, c)),
|
||||||
|
static_cast<T>(values.getValue(1, c)),
|
||||||
|
static_cast<T>(component(pcoords, 0)));
|
||||||
|
auto vbb = internal::lerp(static_cast<T>(values.getValue(3, c)),
|
||||||
|
static_cast<T>(values.getValue(2, c)),
|
||||||
|
static_cast<T>(component(pcoords, 0)));
|
||||||
|
auto vtf = internal::lerp(static_cast<T>(values.getValue(4, c)),
|
||||||
|
static_cast<T>(values.getValue(5, c)),
|
||||||
|
static_cast<T>(component(pcoords, 0)));
|
||||||
|
auto vtb = internal::lerp(static_cast<T>(values.getValue(7, c)),
|
||||||
|
static_cast<T>(values.getValue(6, c)),
|
||||||
|
static_cast<T>(component(pcoords, 0)));
|
||||||
|
auto vb = internal::lerp(vbf, vbb, static_cast<T>(component(pcoords, 1)));
|
||||||
|
auto vt = internal::lerp(vtf, vtb, static_cast<T>(component(pcoords, 1)));
|
||||||
|
auto v = internal::lerp(vb, vt, static_cast<T>(component(pcoords, 2)));
|
||||||
|
component(result, c) = static_cast<ComponentType<Result>>(v);
|
||||||
|
}
|
||||||
|
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
namespace internal
|
||||||
|
{
|
||||||
|
|
||||||
|
template <typename Values, typename CoordType, typename Result>
|
||||||
|
VTKC_EXEC inline void parametricDerivative(Hexahedron,
|
||||||
|
const Values& values,
|
||||||
|
IdComponent comp,
|
||||||
|
const CoordType& pcoords,
|
||||||
|
Result&& result) noexcept
|
||||||
|
{
|
||||||
|
using T = internal::ClosestFloatType<typename Values::ValueType>;
|
||||||
|
T p0 = static_cast<T>(component(pcoords, 0));
|
||||||
|
T p1 = static_cast<T>(component(pcoords, 1));
|
||||||
|
T p2 = static_cast<T>(component(pcoords, 2));
|
||||||
|
T rm = T{1} - p0;
|
||||||
|
T sm = T{1} - p1;
|
||||||
|
T tm = T{1} - p2;
|
||||||
|
|
||||||
|
T dr = (static_cast<T>(values.getValue(0, comp)) * -sm * tm) +
|
||||||
|
(static_cast<T>(values.getValue(1, comp)) * sm * tm) +
|
||||||
|
(static_cast<T>(values.getValue(2, comp)) * p1 * tm) +
|
||||||
|
(static_cast<T>(values.getValue(3, comp)) * -p1 * tm) +
|
||||||
|
(static_cast<T>(values.getValue(4, comp)) * -sm * p2) +
|
||||||
|
(static_cast<T>(values.getValue(5, comp)) * sm * p2) +
|
||||||
|
(static_cast<T>(values.getValue(6, comp)) * p1 * p2) +
|
||||||
|
(static_cast<T>(values.getValue(7, comp)) * -p1 * p2);
|
||||||
|
|
||||||
|
T ds = (static_cast<T>(values.getValue(0, comp)) * -rm * tm) +
|
||||||
|
(static_cast<T>(values.getValue(1, comp)) * -p0 * tm) +
|
||||||
|
(static_cast<T>(values.getValue(2, comp)) * p0 * tm) +
|
||||||
|
(static_cast<T>(values.getValue(3, comp)) * rm * tm) +
|
||||||
|
(static_cast<T>(values.getValue(4, comp)) * -rm * p2) +
|
||||||
|
(static_cast<T>(values.getValue(5, comp)) * -p0 * p2) +
|
||||||
|
(static_cast<T>(values.getValue(6, comp)) * p0 * p2) +
|
||||||
|
(static_cast<T>(values.getValue(7, comp)) * rm * p2);
|
||||||
|
|
||||||
|
T dt = (static_cast<T>(values.getValue(0, comp)) * -rm * sm) +
|
||||||
|
(static_cast<T>(values.getValue(1, comp)) * -p0 * sm) +
|
||||||
|
(static_cast<T>(values.getValue(2, comp)) * -p0 * p1) +
|
||||||
|
(static_cast<T>(values.getValue(3, comp)) * -rm * p1) +
|
||||||
|
(static_cast<T>(values.getValue(4, comp)) * rm * sm) +
|
||||||
|
(static_cast<T>(values.getValue(5, comp)) * p0 * sm) +
|
||||||
|
(static_cast<T>(values.getValue(6, comp)) * p0 * p1) +
|
||||||
|
(static_cast<T>(values.getValue(7, comp)) * rm * p1);
|
||||||
|
|
||||||
|
component(result, 0) = static_cast<ComponentType<Result>>(dr);
|
||||||
|
component(result, 1) = static_cast<ComponentType<Result>>(ds);
|
||||||
|
component(result, 2) = static_cast<ComponentType<Result>>(dt);
|
||||||
|
}
|
||||||
|
|
||||||
|
} // internal
|
||||||
|
|
||||||
|
template <typename Points, typename Values, typename CoordType, typename Result>
|
||||||
|
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<Result>(dx),
|
||||||
|
std::forward<Result>(dy),
|
||||||
|
std::forward<Result>(dz));
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Points, typename PCoordType, typename WCoordType>
|
||||||
|
VTKC_EXEC inline vtkc::ErrorCode parametricToWorld(
|
||||||
|
Hexahedron,
|
||||||
|
const Points& points,
|
||||||
|
const PCoordType& pcoords,
|
||||||
|
WCoordType&& wcoords) noexcept
|
||||||
|
{
|
||||||
|
return interpolate(Hexahedron{}, points, pcoords, std::forward<WCoordType>(wcoords));
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Points, typename WCoordType, typename PCoordType>
|
||||||
|
VTKC_EXEC inline vtkc::ErrorCode worldToParametric(
|
||||||
|
Hexahedron,
|
||||||
|
const Points& points,
|
||||||
|
const WCoordType& wcoords,
|
||||||
|
PCoordType&& pcoords) noexcept
|
||||||
|
{
|
||||||
|
return internal::worldToParametric3D(
|
||||||
|
Hexahedron{}, points, wcoords, std::forward<PCoordType>(pcoords));
|
||||||
|
}
|
||||||
|
|
||||||
|
} // vtkc
|
||||||
|
|
||||||
|
#endif // vtk_c_Hexahedron_h
|
176
vtkc/Line.h
Normal file
176
vtkc/Line.h
Normal file
@ -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 <vtkc/ErrorCode.h>
|
||||||
|
#include <vtkc/Shapes.h>
|
||||||
|
|
||||||
|
#include <vtkc/internal/Common.h>
|
||||||
|
|
||||||
|
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<typename CoordType>
|
||||||
|
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<typename CoordType>
|
||||||
|
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<typename CoordType>
|
||||||
|
VTKC_EXEC inline ComponentType<CoordType> parametricDistance(Line, const CoordType& pcoords) noexcept
|
||||||
|
{
|
||||||
|
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
|
||||||
|
return internal::findParametricDistance(pcoords, 1);
|
||||||
|
}
|
||||||
|
|
||||||
|
template<typename CoordType>
|
||||||
|
VTKC_EXEC inline bool cellInside(Line, const CoordType& pcoords) noexcept
|
||||||
|
{
|
||||||
|
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
|
||||||
|
|
||||||
|
using T = ComponentType<CoordType>;
|
||||||
|
return component(pcoords, 0) >= T{0} && component(pcoords, 0) <= T{1};
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Values, typename CoordType, typename Result>
|
||||||
|
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<typename Values::ValueType>;
|
||||||
|
using ResultCompType = ComponentType<Result>;
|
||||||
|
|
||||||
|
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
|
||||||
|
{
|
||||||
|
auto ival = internal::lerp(static_cast<ProcessingType>(values.getValue(0, c)),
|
||||||
|
static_cast<ProcessingType>(values.getValue(1, c)),
|
||||||
|
static_cast<ProcessingType>(component(pcoords, 0)));
|
||||||
|
component(result, c) = static_cast<ResultCompType>(ival);
|
||||||
|
}
|
||||||
|
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Points, typename Values, typename CoordType, typename Result>
|
||||||
|
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<typename Values::ValueType>;
|
||||||
|
using ResultCompType = ComponentType<Result>;
|
||||||
|
|
||||||
|
ProcessingType dPt[3] = {
|
||||||
|
static_cast<ProcessingType>(points.getValue(1, 0) - points.getValue(0, 0)),
|
||||||
|
static_cast<ProcessingType>(points.getValue(1, 1) - points.getValue(0, 1)),
|
||||||
|
static_cast<ProcessingType>(points.getValue(1, 2) - points.getValue(0, 2)) };
|
||||||
|
|
||||||
|
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
|
||||||
|
{
|
||||||
|
auto dv = static_cast<ProcessingType>(values.getValue(1, c) - values.getValue(0, c));
|
||||||
|
component(dx, c) = (dPt[0] != 0.0f) ? static_cast<ResultCompType>(dv/dPt[0]) : ResultCompType{0};
|
||||||
|
component(dy, c) = (dPt[1] != 0.0f) ? static_cast<ResultCompType>(dv/dPt[1]) : ResultCompType{0};
|
||||||
|
component(dz, c) = (dPt[2] != 0.0f) ? static_cast<ResultCompType>(dv/dPt[2]) : ResultCompType{0};
|
||||||
|
}
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Points, typename PCoordType, typename WCoordType>
|
||||||
|
VTKC_EXEC inline vtkc::ErrorCode parametricToWorld(
|
||||||
|
Line,
|
||||||
|
const Points& points,
|
||||||
|
const PCoordType& pcoords,
|
||||||
|
WCoordType&& wcoords) noexcept
|
||||||
|
{
|
||||||
|
return interpolate(Line{}, points, pcoords, std::forward<WCoordType>(wcoords));
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Points, typename WCoordType, typename PCoordType>
|
||||||
|
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<PCoordType>;
|
||||||
|
internal::Vector<T, 3> p0(static_cast<T>(points.getValue(0, 0)),
|
||||||
|
static_cast<T>(points.getValue(0, 1)),
|
||||||
|
static_cast<T>(points.getValue(0, 2)));
|
||||||
|
internal::Vector<T, 3> p1(static_cast<T>(points.getValue(1, 0)),
|
||||||
|
static_cast<T>(points.getValue(1, 1)),
|
||||||
|
static_cast<T>(points.getValue(1, 2)));
|
||||||
|
internal::Vector<T, 3> wc(static_cast<T>(component(wcoords, 0)),
|
||||||
|
static_cast<T>(component(wcoords, 1)),
|
||||||
|
static_cast<T>(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
|
178
vtkc/Pixel.h
Normal file
178
vtkc/Pixel.h
Normal file
@ -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 <vtkc/ErrorCode.h>
|
||||||
|
#include <vtkc/Quad.h>
|
||||||
|
|
||||||
|
#include <vtkc/internal/Common.h>
|
||||||
|
|
||||||
|
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 <typename Points, typename T>
|
||||||
|
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<T>(points.getValue(2, i) - points.getValue(0, i));
|
||||||
|
if (spacing[i] == T{0})
|
||||||
|
{
|
||||||
|
zeros |= 1 << i;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return zeros;
|
||||||
|
}
|
||||||
|
|
||||||
|
} // internal
|
||||||
|
|
||||||
|
template <typename Points, typename Values, typename CoordType, typename Result>
|
||||||
|
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<typename Values::ValueType>;
|
||||||
|
using ResultCompType = ComponentType<Result>;
|
||||||
|
|
||||||
|
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<ResultCompType>(dvdp[0] / spacing[1]);
|
||||||
|
component(dz, c) = static_cast<ResultCompType>(dvdp[1] / spacing[2]);
|
||||||
|
break;
|
||||||
|
case 2: // xz plane
|
||||||
|
component(dx, c) = static_cast<ResultCompType>(dvdp[0] / spacing[0]);
|
||||||
|
component(dy, c) = ResultCompType{0};
|
||||||
|
component(dz, c) = static_cast<ResultCompType>(dvdp[1] / spacing[2]);
|
||||||
|
break;
|
||||||
|
case 4: // xy plane
|
||||||
|
component(dx, c) = static_cast<ResultCompType>(dvdp[0] / spacing[0]);
|
||||||
|
component(dy, c) = static_cast<ResultCompType>(dvdp[1] / spacing[1]);
|
||||||
|
component(dz, c) = ResultCompType{0};
|
||||||
|
break;
|
||||||
|
default:
|
||||||
|
return ErrorCode::DEGENERATE_CELL_DETECTED;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Points, typename PCoordType, typename WCoordType>
|
||||||
|
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<T>(component(pcoords, 0)));
|
||||||
|
component(wcoords, 2) = points.getValue(0, 2) +
|
||||||
|
(spacing[2] * static_cast<T>(component(pcoords, 1)));
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
case 2: // xz plane
|
||||||
|
component(wcoords, 0) = points.getValue(0, 0) +
|
||||||
|
(spacing[0] * static_cast<T>(component(pcoords, 0)));
|
||||||
|
component(wcoords, 1) = points.getValue(0, 1);
|
||||||
|
component(wcoords, 2) = points.getValue(0, 2) +
|
||||||
|
(spacing[2] * static_cast<T>(component(pcoords, 1)));
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
case 4: // xy plane
|
||||||
|
component(wcoords, 0) = points.getValue(0, 0) +
|
||||||
|
(spacing[0] * static_cast<T>(component(pcoords, 0)));
|
||||||
|
component(wcoords, 1) = points.getValue(0, 1) +
|
||||||
|
(spacing[1] * static_cast<T>(component(pcoords, 1)));
|
||||||
|
component(wcoords, 2) = points.getValue(0, 2);
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
default:
|
||||||
|
return ErrorCode::DEGENERATE_CELL_DETECTED;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Points, typename WCoordType, typename PCoordType>
|
||||||
|
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<PCoordType>;
|
||||||
|
|
||||||
|
T spacing[3];
|
||||||
|
int zeros = internal::getPixelSpacing(points, spacing);
|
||||||
|
|
||||||
|
switch (zeros)
|
||||||
|
{
|
||||||
|
case 1: // yz plane
|
||||||
|
component(pcoords, 0) = static_cast<T>(component(wcoords, 1) - points.getValue(0, 1)) /
|
||||||
|
spacing[1];
|
||||||
|
component(pcoords, 1) = static_cast<T>(component(wcoords, 2) - points.getValue(0, 2)) /
|
||||||
|
spacing[2];
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
case 2: // xz plane
|
||||||
|
component(pcoords, 0) = static_cast<T>(component(wcoords, 0) - points.getValue(0, 0)) /
|
||||||
|
spacing[0];
|
||||||
|
component(pcoords, 1) = static_cast<T>(component(wcoords, 2) - points.getValue(0, 2)) /
|
||||||
|
spacing[2];
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
case 4: // xy plane
|
||||||
|
component(pcoords, 0) = static_cast<T>(component(wcoords, 0) - points.getValue(0, 0)) /
|
||||||
|
spacing[0];
|
||||||
|
component(pcoords, 1) = static_cast<T>(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
|
541
vtkc/Polygon.h
Normal file
541
vtkc/Polygon.h
Normal file
@ -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 <vtkc/ErrorCode.h>
|
||||||
|
#include <vtkc/Quad.h>
|
||||||
|
#include <vtkc/Shapes.h>
|
||||||
|
#include <vtkc/Triangle.h>
|
||||||
|
|
||||||
|
#include <vtkc/internal/Common.h>
|
||||||
|
|
||||||
|
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<typename CoordType>
|
||||||
|
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<typename CoordType>
|
||||||
|
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<CoordType>;
|
||||||
|
constexpr double two_pi = 2.0 * 3.14159265359;
|
||||||
|
auto angle = (static_cast<T>(pointId) * static_cast<T>(two_pi)) / static_cast<T>(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<typename CoordType>
|
||||||
|
VTKC_EXEC inline ComponentType<CoordType> 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<typename CoordType>
|
||||||
|
VTKC_EXEC inline bool cellInside(Polygon tag, const CoordType& pcoords) noexcept
|
||||||
|
{
|
||||||
|
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
|
||||||
|
|
||||||
|
using T = ComponentType<CoordType>;
|
||||||
|
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 <typename CoordType>
|
||||||
|
VTKC_EXEC inline vtkc::ErrorCode polygonToSubTrianglePCoords(
|
||||||
|
Polygon tag,
|
||||||
|
const CoordType& polygonPC,
|
||||||
|
IdComponent& p0,
|
||||||
|
IdComponent& p1,
|
||||||
|
ComponentType<CoordType> trianglePC[2]) noexcept
|
||||||
|
{
|
||||||
|
using T = ComponentType<CoordType>;
|
||||||
|
|
||||||
|
constexpr T epsilon = std::is_same<T, float>::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<T>(two_pi);
|
||||||
|
}
|
||||||
|
T deltaAngle = static_cast<T>(two_pi) / static_cast<T>(tag.numberOfPoints());
|
||||||
|
|
||||||
|
p0 = static_cast<IdComponent>(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 <typename Values>
|
||||||
|
VTKC_EXEC inline typename Values::ValueType polygonInterpolateComponentAtCenter(
|
||||||
|
Polygon tag, const Values& values, IdComponent comp) noexcept
|
||||||
|
{
|
||||||
|
using T = internal::ClosestFloatType<typename Values::ValueType>;
|
||||||
|
|
||||||
|
auto weight = T{1} / static_cast<T>(tag.numberOfPoints());
|
||||||
|
auto result = static_cast<T>(values.getValue(0, comp));
|
||||||
|
for (IdComponent i = 1; i < tag.numberOfPoints(); ++i)
|
||||||
|
{
|
||||||
|
result += static_cast<T>(values.getValue(i, comp));
|
||||||
|
}
|
||||||
|
result *= weight;
|
||||||
|
|
||||||
|
return static_cast<typename Values::ValueType>(result);
|
||||||
|
}
|
||||||
|
|
||||||
|
} // namespace internal
|
||||||
|
|
||||||
|
template <typename Values, typename CoordType, typename Result>
|
||||||
|
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>(result));
|
||||||
|
case 4:
|
||||||
|
return interpolate(Quad{}, values, pcoords, std::forward<Result>(result));
|
||||||
|
default:
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
|
||||||
|
using ResultCompType = ComponentType<Result>;
|
||||||
|
using ProcessingType = internal::ClosestFloatType<typename Values::ValueType>;
|
||||||
|
|
||||||
|
IdComponent p0, p1;
|
||||||
|
ComponentType<CoordType> 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<ProcessingType>(internal::polygonInterpolateComponentAtCenter(tag, values, c));
|
||||||
|
triVals[1] = static_cast<ProcessingType>(values.getValue(p0, c));
|
||||||
|
triVals[2] = static_cast<ProcessingType>(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 <typename CoordType>
|
||||||
|
VTKC_EXEC inline void polygonGetTriangleAroundPCoords(
|
||||||
|
const CoordType& pcoords, ComponentType<CoordType> pc1[2], ComponentType<CoordType> pc2[2]) noexcept
|
||||||
|
{
|
||||||
|
using T = ComponentType<CoordType>;
|
||||||
|
|
||||||
|
// Find the unit vector pointing from the center of the polygon to pcoords
|
||||||
|
Vector<T, 2> 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 <typename Points, typename Values, typename CoordType, typename Result>
|
||||||
|
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<Result>(dx),
|
||||||
|
std::forward<Result>(dy),
|
||||||
|
std::forward<Result>(dz));
|
||||||
|
case 4:
|
||||||
|
return derivative(Quad{},
|
||||||
|
points,
|
||||||
|
values,
|
||||||
|
pcoords,
|
||||||
|
std::forward<Result>(dx),
|
||||||
|
std::forward<Result>(dy),
|
||||||
|
std::forward<Result>(dz));
|
||||||
|
default:
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
|
||||||
|
using ResultCompType = ComponentType<Result>;
|
||||||
|
using ProcessingType = internal::ClosestFloatType<typename Values::ValueType>;
|
||||||
|
|
||||||
|
// Get the parametric coordinates of a small triangle, with pcoords as one of the vertices
|
||||||
|
ComponentType<CoordType> ptPc1[2], ptPc2[2];
|
||||||
|
internal::polygonGetTriangleAroundPCoords(pcoords, ptPc1, ptPc2);
|
||||||
|
|
||||||
|
// Compute world coordinates of the points of the triangle
|
||||||
|
internal::Vector<ProcessingType, 3> 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<ProcessingType> triSpace(triPts[0], triPts[1], triPts[2]);
|
||||||
|
internal::Vector<ProcessingType, 2> pts2d[3];
|
||||||
|
for (int i = 0; i < 3; ++i)
|
||||||
|
{
|
||||||
|
pts2d[i] = triSpace.to2DPoint(triPts[i]);
|
||||||
|
}
|
||||||
|
|
||||||
|
// pre-compute once
|
||||||
|
internal::Matrix<ProcessingType, 2, 2> jacobian;
|
||||||
|
internal::jacobian2D(Triangle{}, makeFieldAccessorNestedSOA(pts2d, 2), nullptr, jacobian);
|
||||||
|
internal::Matrix<ProcessingType, 2, 2> 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<CoordType> 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<ProcessingType>(internal::polygonInterpolateComponentAtCenter(tag, values, c));
|
||||||
|
ProcessingType triVals[3];
|
||||||
|
for (int i = 0; i < 3; ++i)
|
||||||
|
{
|
||||||
|
ProcessingType field[3] = {vCenter,
|
||||||
|
static_cast<ProcessingType>(values.getValue(subP1P2[i][0], c)),
|
||||||
|
static_cast<ProcessingType>(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<ProcessingType, 2> dvdp;
|
||||||
|
parametricDerivative(Triangle{}, makeFieldAccessorNestedSOA(triVals), 0, nullptr, dvdp);
|
||||||
|
auto d2D = matrixMultiply(dvdp, invJacobian);
|
||||||
|
auto d3D = triSpace.to3DVec(d2D);
|
||||||
|
|
||||||
|
component(dx, c) = static_cast<ResultCompType>(d3D[0]);
|
||||||
|
component(dy, c) = static_cast<ResultCompType>(d3D[1]);
|
||||||
|
component(dz, c) = static_cast<ResultCompType>(d3D[2]);
|
||||||
|
}
|
||||||
|
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Points, typename PCoordType, typename WCoordType>
|
||||||
|
VTKC_EXEC inline vtkc::ErrorCode parametricToWorld(
|
||||||
|
Polygon tag,
|
||||||
|
const Points& points,
|
||||||
|
const PCoordType& pcoords,
|
||||||
|
WCoordType&& wcoords) noexcept
|
||||||
|
{
|
||||||
|
return interpolate(tag, points, pcoords, std::forward<WCoordType>(wcoords));
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Points, typename WCoordType, typename PCoordType>
|
||||||
|
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<WCoordType>;
|
||||||
|
auto numPoints = tag.numberOfPoints();
|
||||||
|
|
||||||
|
// Find the position of the center point.
|
||||||
|
internal::Vector<T, 3> 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<T>(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<T, 3> polygonNormal;
|
||||||
|
{
|
||||||
|
internal::Vector<T, 3> v1p1, v1p2;
|
||||||
|
points.getTuple(0, v1p1);
|
||||||
|
points.getTuple(numPoints / 3, v1p2);
|
||||||
|
|
||||||
|
internal::Vector<T, 3> 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<T, 3> wc{ component(wcoords, 0), component(wcoords, 1), component(wcoords, 2) };
|
||||||
|
IdComponent firstPointIndex;
|
||||||
|
IdComponent secondPointIndex = 0;
|
||||||
|
internal::Vector<T, 3> 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<T, 3> triangleWCoords[3] = { wcoordCenter, firstPoint, secondPoint };
|
||||||
|
internal::Vector<T, 3> 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<ComponentType<PCoordType>>(wc[0]);
|
||||||
|
component(pcoords, 1) = static_cast<ComponentType<PCoordType>>(wc[1]);
|
||||||
|
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
} //namespace vtkc
|
||||||
|
|
||||||
|
#endif //vtk_c_Polygon_h
|
293
vtkc/Pyramid.h
Normal file
293
vtkc/Pyramid.h
Normal file
@ -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 <vtkc/ErrorCode.h>
|
||||||
|
#include <vtkc/Shapes.h>
|
||||||
|
|
||||||
|
#include <vtkc/internal/Common.h>
|
||||||
|
|
||||||
|
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<typename CoordType>
|
||||||
|
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<typename CoordType>
|
||||||
|
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<typename CoordType>
|
||||||
|
VTKC_EXEC inline ComponentType<CoordType> parametricDistance(Pyramid, const CoordType& pcoords) noexcept
|
||||||
|
{
|
||||||
|
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
|
||||||
|
return internal::findParametricDistance(pcoords, 3);
|
||||||
|
}
|
||||||
|
|
||||||
|
template<typename CoordType>
|
||||||
|
VTKC_EXEC inline bool cellInside(Pyramid, const CoordType& pcoords) noexcept
|
||||||
|
{
|
||||||
|
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
|
||||||
|
|
||||||
|
using T = ComponentType<CoordType>;
|
||||||
|
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 <typename Values, typename CoordType, typename Result>
|
||||||
|
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<typename Values::ValueType>;
|
||||||
|
|
||||||
|
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
|
||||||
|
{
|
||||||
|
auto baseV0 = internal::lerp(static_cast<T>(values.getValue(0, c)),
|
||||||
|
static_cast<T>(values.getValue(1, c)),
|
||||||
|
static_cast<T>(component(pcoords, 0)));
|
||||||
|
auto baseV1 = internal::lerp(static_cast<T>(values.getValue(3, c)),
|
||||||
|
static_cast<T>(values.getValue(2, c)),
|
||||||
|
static_cast<T>(component(pcoords, 0)));
|
||||||
|
|
||||||
|
auto baseV = internal::lerp(baseV0, baseV1, static_cast<T>(component(pcoords, 1)));
|
||||||
|
|
||||||
|
auto v = internal::lerp(baseV,
|
||||||
|
static_cast<T>(values.getValue(4, c)),
|
||||||
|
static_cast<T>(component(pcoords, 2)));
|
||||||
|
|
||||||
|
component(result, c) = static_cast<ComponentType<Result>>(v);
|
||||||
|
}
|
||||||
|
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
namespace internal
|
||||||
|
{
|
||||||
|
|
||||||
|
template <typename Values, typename CoordType, typename Result>
|
||||||
|
VTKC_EXEC inline void parametricDerivative(Pyramid,
|
||||||
|
const Values& values,
|
||||||
|
IdComponent comp,
|
||||||
|
const CoordType& pcoords,
|
||||||
|
Result&& result) noexcept
|
||||||
|
{
|
||||||
|
using T = internal::ClosestFloatType<typename Values::ValueType>;
|
||||||
|
T p0 = static_cast<T>(component(pcoords, 0));
|
||||||
|
T p1 = static_cast<T>(component(pcoords, 1));
|
||||||
|
T p2 = static_cast<T>(component(pcoords, 2));
|
||||||
|
T rm = T{1} - p0;
|
||||||
|
T sm = T{1} - p1;
|
||||||
|
T tm = T{1} - p2;
|
||||||
|
|
||||||
|
T dr = (static_cast<T>(values.getValue(0, comp)) * -sm * tm) +
|
||||||
|
(static_cast<T>(values.getValue(1, comp)) * sm * tm) +
|
||||||
|
(static_cast<T>(values.getValue(2, comp)) * p1 * tm) +
|
||||||
|
(static_cast<T>(values.getValue(3, comp)) * -p1 * tm);
|
||||||
|
|
||||||
|
T ds = (static_cast<T>(values.getValue(0, comp)) * -rm * tm) +
|
||||||
|
(static_cast<T>(values.getValue(1, comp)) * -p0 * tm) +
|
||||||
|
(static_cast<T>(values.getValue(2, comp)) * p0 * tm) +
|
||||||
|
(static_cast<T>(values.getValue(3, comp)) * rm * tm);
|
||||||
|
|
||||||
|
T dt = (static_cast<T>(values.getValue(0, comp)) * -rm * sm) +
|
||||||
|
(static_cast<T>(values.getValue(1, comp)) * -p0 * sm) +
|
||||||
|
(static_cast<T>(values.getValue(2, comp)) * -p0 * p1) +
|
||||||
|
(static_cast<T>(values.getValue(3, comp)) * -rm * p1) +
|
||||||
|
(static_cast<T>(values.getValue(4, comp)));
|
||||||
|
|
||||||
|
component(result, 0) = static_cast<ComponentType<Result>>(dr);
|
||||||
|
component(result, 1) = static_cast<ComponentType<Result>>(ds);
|
||||||
|
component(result, 2) = static_cast<ComponentType<Result>>(dt);
|
||||||
|
}
|
||||||
|
|
||||||
|
} // internal
|
||||||
|
|
||||||
|
template <typename Points, typename Values, typename CoordType, typename Result>
|
||||||
|
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<typename Values::ValueType>;
|
||||||
|
using ResultCompType = ComponentType<Result>;
|
||||||
|
|
||||||
|
if (component(pcoords, 2) > ComponentType<CoordType>(.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<ProcessingType, 3, 3> j, ij1, ij2;
|
||||||
|
|
||||||
|
const ComponentType<CoordType> 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<CoordType> 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<ProcessingType, 3> 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<ResultCompType>((d2[0] * 2.0f) - d1[0]);
|
||||||
|
component(dy, c) = static_cast<ResultCompType>((d2[1] * 2.0f) - d1[1]);
|
||||||
|
component(dz, c) = static_cast<ResultCompType>((d2[2] * 2.0f) - d1[2]);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
return internal::derivative3D(Pyramid{},
|
||||||
|
points,
|
||||||
|
values,
|
||||||
|
pcoords,
|
||||||
|
std::forward<Result>(dx),
|
||||||
|
std::forward<Result>(dy),
|
||||||
|
std::forward<Result>(dz));
|
||||||
|
}
|
||||||
|
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Points, typename PCoordType, typename WCoordType>
|
||||||
|
VTKC_EXEC inline vtkc::ErrorCode parametricToWorld(
|
||||||
|
Pyramid,
|
||||||
|
const Points& points,
|
||||||
|
const PCoordType& pcoords,
|
||||||
|
WCoordType&& wcoords) noexcept
|
||||||
|
{
|
||||||
|
return interpolate(Pyramid{}, points, pcoords, std::forward<WCoordType>(wcoords));
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Points, typename WCoordType, typename PCoordType>
|
||||||
|
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<PCoordType>;
|
||||||
|
|
||||||
|
internal::Vector<TIn, 3> 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<TOut, 3> pcBaseCenter(0.5f, 0.5f, 0.0f);
|
||||||
|
internal::Vector<TIn, 3> 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<PCoordType>(pcoords));
|
||||||
|
}
|
||||||
|
|
||||||
|
} // vtkc
|
||||||
|
|
||||||
|
#endif // vtk_c_Pyramid_h
|
195
vtkc/Quad.h
Normal file
195
vtkc/Quad.h
Normal file
@ -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 <vtkc/ErrorCode.h>
|
||||||
|
#include <vtkc/Shapes.h>
|
||||||
|
|
||||||
|
#include <vtkc/internal/Common.h>
|
||||||
|
|
||||||
|
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<typename CoordType>
|
||||||
|
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<typename CoordType>
|
||||||
|
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<typename CoordType>
|
||||||
|
VTKC_EXEC inline ComponentType<CoordType> parametricDistance(Quad, const CoordType& pcoords) noexcept
|
||||||
|
{
|
||||||
|
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
|
||||||
|
return internal::findParametricDistance(pcoords, 2);
|
||||||
|
}
|
||||||
|
|
||||||
|
template<typename CoordType>
|
||||||
|
VTKC_EXEC inline bool cellInside(Quad, const CoordType& pcoords) noexcept
|
||||||
|
{
|
||||||
|
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
|
||||||
|
|
||||||
|
using T = ComponentType<CoordType>;
|
||||||
|
return component(pcoords, 0) >= T{0} && component(pcoords, 0) <= T{1} &&
|
||||||
|
component(pcoords, 1) >= T{0} && component(pcoords, 1) <= T{1};
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Values, typename CoordType, typename Result>
|
||||||
|
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<typename Values::ValueType>;
|
||||||
|
|
||||||
|
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
|
||||||
|
{
|
||||||
|
auto v0 = internal::lerp(static_cast<T>(values.getValue(0, c)),
|
||||||
|
static_cast<T>(values.getValue(1, c)),
|
||||||
|
static_cast<T>(component(pcoords, 0)));
|
||||||
|
auto v1 = internal::lerp(static_cast<T>(values.getValue(3, c)),
|
||||||
|
static_cast<T>(values.getValue(2, c)),
|
||||||
|
static_cast<T>(component(pcoords, 0)));
|
||||||
|
auto v = internal::lerp(v0, v1, static_cast<T>(component(pcoords, 1)));
|
||||||
|
component(result, c) = static_cast<ComponentType<Result>>(v);
|
||||||
|
}
|
||||||
|
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
namespace internal
|
||||||
|
{
|
||||||
|
|
||||||
|
template <typename Values, typename CoordType, typename Result>
|
||||||
|
VTKC_EXEC inline void parametricDerivative(
|
||||||
|
Quad, const Values& values, IdComponent comp, const CoordType& pcoords, Result&& result) noexcept
|
||||||
|
{
|
||||||
|
using T = internal::ClosestFloatType<typename Values::ValueType>;
|
||||||
|
T p0 = static_cast<T>(component(pcoords, 0));
|
||||||
|
T p1 = static_cast<T>(component(pcoords, 1));
|
||||||
|
T rm = T{1} - p0;
|
||||||
|
T sm = T{1} - p1;
|
||||||
|
|
||||||
|
T dr = (static_cast<T>(values.getValue(0, comp)) * -sm) +
|
||||||
|
(static_cast<T>(values.getValue(1, comp)) * sm) +
|
||||||
|
(static_cast<T>(values.getValue(2, comp)) * p1) +
|
||||||
|
(static_cast<T>(values.getValue(3, comp)) * -p1);
|
||||||
|
T ds = (static_cast<T>(values.getValue(0, comp)) * -rm) +
|
||||||
|
(static_cast<T>(values.getValue(1, comp)) * -p0) +
|
||||||
|
(static_cast<T>(values.getValue(2, comp)) * p0) +
|
||||||
|
(static_cast<T>(values.getValue(3, comp)) * rm);
|
||||||
|
|
||||||
|
component(result, 0) = static_cast<ComponentType<Result>>(dr);
|
||||||
|
component(result, 1) = static_cast<ComponentType<Result>>(ds);
|
||||||
|
}
|
||||||
|
|
||||||
|
} // internal
|
||||||
|
|
||||||
|
template <typename Points, typename Values, typename CoordType, typename Result>
|
||||||
|
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<Result>(dx),
|
||||||
|
std::forward<Result>(dy),
|
||||||
|
std::forward<Result>(dz));
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Points, typename PCoordType, typename WCoordType>
|
||||||
|
VTKC_EXEC inline vtkc::ErrorCode parametricToWorld(
|
||||||
|
Quad,
|
||||||
|
const Points& points,
|
||||||
|
const PCoordType& pcoords,
|
||||||
|
WCoordType&& wcoords) noexcept
|
||||||
|
{
|
||||||
|
return interpolate(Quad{}, points, pcoords, std::forward<WCoordType>(wcoords));
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Points, typename WCoordType, typename PCoordType>
|
||||||
|
VTKC_EXEC inline vtkc::ErrorCode worldToParametric(
|
||||||
|
Quad,
|
||||||
|
const Points& points,
|
||||||
|
const WCoordType& wcoords,
|
||||||
|
PCoordType&& pcoords) noexcept
|
||||||
|
{
|
||||||
|
return internal::worldToParametric2D(Quad{}, points, wcoords, std::forward<PCoordType>(pcoords));
|
||||||
|
}
|
||||||
|
|
||||||
|
} //namespace vtkc
|
||||||
|
|
||||||
|
#endif //vtk_c_Quad_h
|
142
vtkc/Shapes.h
Normal file
142
vtkc/Shapes.h
Normal file
@ -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 <vtkc/internal/Config.h>
|
||||||
|
|
||||||
|
#include <cstdint>
|
||||||
|
|
||||||
|
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
|
208
vtkc/Tetra.h
Normal file
208
vtkc/Tetra.h
Normal file
@ -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 <vtkc/ErrorCode.h>
|
||||||
|
#include <vtkc/Shapes.h>
|
||||||
|
|
||||||
|
#include <vtkc/internal/Common.h>
|
||||||
|
|
||||||
|
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<typename CoordType>
|
||||||
|
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<typename CoordType>
|
||||||
|
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<typename CoordType>
|
||||||
|
VTKC_EXEC inline ComponentType<CoordType> parametricDistance(Tetra, const CoordType& pcoords) noexcept
|
||||||
|
{
|
||||||
|
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
|
||||||
|
|
||||||
|
ComponentType<CoordType> weights[4];
|
||||||
|
weights[0] = ComponentType<CoordType>{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<typename CoordType>
|
||||||
|
VTKC_EXEC inline bool cellInside(Tetra, const CoordType& pcoords) noexcept
|
||||||
|
{
|
||||||
|
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
|
||||||
|
|
||||||
|
using T = ComponentType<CoordType>;
|
||||||
|
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 <typename CoordType, typename T>
|
||||||
|
VTKC_EXEC inline void interpolationWeights(Tetra, const CoordType& pcoords, T weights[4]) noexcept
|
||||||
|
{
|
||||||
|
weights[1] = static_cast<T>(component(pcoords, 0));
|
||||||
|
weights[2] = static_cast<T>(component(pcoords, 1));
|
||||||
|
weights[3] = static_cast<T>(component(pcoords, 2));
|
||||||
|
weights[0] = T(1) - weights[1] - weights[2] - weights[3];
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Values, typename CoordType, typename Result>
|
||||||
|
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<typename Values::ValueType>;
|
||||||
|
|
||||||
|
auto w0 = T(1) - static_cast<T>(component(pcoords, 0) + component(pcoords, 1) + component(pcoords, 2));
|
||||||
|
auto w1 = static_cast<T>(component(pcoords, 0));
|
||||||
|
auto w2 = static_cast<T>(component(pcoords, 1));
|
||||||
|
auto w3 = static_cast<T>(component(pcoords, 2));
|
||||||
|
|
||||||
|
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
|
||||||
|
{
|
||||||
|
auto v = static_cast<T>(values.getValue(0, c)) * w0 +
|
||||||
|
static_cast<T>(values.getValue(1, c)) * w1 +
|
||||||
|
static_cast<T>(values.getValue(2, c)) * w2 +
|
||||||
|
static_cast<T>(values.getValue(3, c)) * w3;
|
||||||
|
component(result, c) = static_cast<ComponentType<Result>>(v);
|
||||||
|
}
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
namespace internal
|
||||||
|
{
|
||||||
|
|
||||||
|
template <typename Values, typename CoordType, typename Result>
|
||||||
|
VTKC_EXEC inline void parametricDerivative(
|
||||||
|
Tetra, const Values& values, IdComponent comp, const CoordType&, Result&& result) noexcept
|
||||||
|
{
|
||||||
|
component(result, 0) = static_cast<ComponentType<Result>>(values.getValue(1, comp) -
|
||||||
|
values.getValue(0, comp));
|
||||||
|
component(result, 1) = static_cast<ComponentType<Result>>(values.getValue(2, comp) -
|
||||||
|
values.getValue(0, comp));
|
||||||
|
component(result, 2) = static_cast<ComponentType<Result>>(values.getValue(3, comp) -
|
||||||
|
values.getValue(0, comp));
|
||||||
|
}
|
||||||
|
|
||||||
|
} // internal
|
||||||
|
|
||||||
|
template <typename Points, typename Values, typename CoordType, typename Result>
|
||||||
|
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<Result>(dx),
|
||||||
|
std::forward<Result>(dy),
|
||||||
|
std::forward<Result>(dz));
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Points, typename PCoordType, typename WCoordType>
|
||||||
|
VTKC_EXEC inline vtkc::ErrorCode parametricToWorld(
|
||||||
|
Tetra,
|
||||||
|
const Points& points,
|
||||||
|
const PCoordType& pcoords,
|
||||||
|
WCoordType&& wcoords) noexcept
|
||||||
|
{
|
||||||
|
return interpolate(Tetra{}, points, pcoords, std::forward<WCoordType>(wcoords));
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Points, typename WCoordType, typename PCoordType>
|
||||||
|
VTKC_EXEC inline vtkc::ErrorCode worldToParametric(
|
||||||
|
Tetra,
|
||||||
|
const Points& points,
|
||||||
|
const WCoordType& wcoords,
|
||||||
|
PCoordType&& pcoords) noexcept
|
||||||
|
{
|
||||||
|
return internal::worldToParametric3D(Tetra{}, points, wcoords, std::forward<PCoordType>(pcoords));
|
||||||
|
}
|
||||||
|
|
||||||
|
} // vtkc
|
||||||
|
|
||||||
|
#endif // vtk_c_Tetra_h
|
271
vtkc/Triangle.h
Normal file
271
vtkc/Triangle.h
Normal file
@ -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 <vtkc/ErrorCode.h>
|
||||||
|
#include <vtkc/Shapes.h>
|
||||||
|
|
||||||
|
#include <vtkc/internal/Common.h>
|
||||||
|
|
||||||
|
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<typename CoordType>
|
||||||
|
VTKC_EXEC inline vtkc::ErrorCode parametricCenter(Triangle, CoordType&& pcoords) noexcept
|
||||||
|
{
|
||||||
|
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
|
||||||
|
|
||||||
|
using T = ComponentType<CoordType>;
|
||||||
|
component(pcoords, 0) = T(1)/T(3);
|
||||||
|
component(pcoords, 1) = T(1)/T(3);
|
||||||
|
component(pcoords, 2) = T(0);
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
template<typename CoordType>
|
||||||
|
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<typename CoordType>
|
||||||
|
VTKC_EXEC inline ComponentType<CoordType> parametricDistance(Triangle, const CoordType& pcoords) noexcept
|
||||||
|
{
|
||||||
|
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
|
||||||
|
|
||||||
|
ComponentType<CoordType> weights[3];
|
||||||
|
weights[0] = ComponentType<CoordType>{1} - component(pcoords, 0) - component(pcoords, 1);
|
||||||
|
weights[1] = component(pcoords, 0);
|
||||||
|
weights[2] = component(pcoords, 1);
|
||||||
|
return internal::findParametricDistance(weights, 3);
|
||||||
|
}
|
||||||
|
|
||||||
|
template<typename CoordType>
|
||||||
|
VTKC_EXEC inline bool cellInside(Triangle, const CoordType& pcoords) noexcept
|
||||||
|
{
|
||||||
|
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
|
||||||
|
|
||||||
|
using T = ComponentType<CoordType>;
|
||||||
|
return component(pcoords, 0) >= T{0} && component(pcoords, 1) >= T{0} &&
|
||||||
|
(component(pcoords, 0) + component(pcoords, 1)) <= T{1};
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Values, typename CoordType, typename Result>
|
||||||
|
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<typename Values::ValueType>;
|
||||||
|
|
||||||
|
auto w0 = T(1) - static_cast<T>(component(pcoords, 0) + component(pcoords, 1));
|
||||||
|
auto w1 = static_cast<T>(component(pcoords, 0));
|
||||||
|
auto w2 = static_cast<T>(component(pcoords, 1));
|
||||||
|
|
||||||
|
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
|
||||||
|
{
|
||||||
|
auto v = static_cast<T>(values.getValue(0, c)) * w0 +
|
||||||
|
static_cast<T>(values.getValue(1, c)) * w1 +
|
||||||
|
static_cast<T>(values.getValue(2, c)) * w2;
|
||||||
|
component(result, c) = static_cast<ComponentType<Result>>(v);
|
||||||
|
}
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
namespace internal
|
||||||
|
{
|
||||||
|
|
||||||
|
template <typename Values, typename CoordType, typename Result>
|
||||||
|
VTKC_EXEC inline void parametricDerivative(
|
||||||
|
Triangle, const Values& values, IdComponent comp, const CoordType&, Result&& result) noexcept
|
||||||
|
{
|
||||||
|
component(result, 0) = static_cast<ComponentType<Result>>(values.getValue(1, comp) -
|
||||||
|
values.getValue(0, comp));
|
||||||
|
component(result, 1) = static_cast<ComponentType<Result>>(values.getValue(2, comp) -
|
||||||
|
values.getValue(0, comp));
|
||||||
|
}
|
||||||
|
|
||||||
|
} // internal
|
||||||
|
|
||||||
|
template <typename Points, typename Values, typename CoordType, typename Result>
|
||||||
|
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<Result>(dx),
|
||||||
|
std::forward<Result>(dy),
|
||||||
|
std::forward<Result>(dz));
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Points, typename PCoordType, typename WCoordType>
|
||||||
|
VTKC_EXEC inline vtkc::ErrorCode parametricToWorld(
|
||||||
|
Triangle,
|
||||||
|
const Points& points,
|
||||||
|
const PCoordType& pcoords,
|
||||||
|
WCoordType&& wcoords) noexcept
|
||||||
|
{
|
||||||
|
return interpolate(Triangle{}, points, pcoords, std::forward<WCoordType>(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 <typename Points, typename WCoordType, typename PCoordType>
|
||||||
|
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<PCoordType>;
|
||||||
|
|
||||||
|
internal::Vector<TIn, 3> pts[3];
|
||||||
|
for (int i = 0; i < 3; ++i)
|
||||||
|
{
|
||||||
|
points.getTuple(i, pts[i]);
|
||||||
|
}
|
||||||
|
|
||||||
|
internal::Vector<TIn, 3> 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<TOut>(internal::dot(wc - p0, planeNormal)) /
|
||||||
|
static_cast<TOut>(internal::dot(p1 - p0, planeNormal));
|
||||||
|
}
|
||||||
|
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
} //namespace vtkc
|
||||||
|
|
||||||
|
#endif //vtk_c_Triangle_h
|
131
vtkc/Vertex.h
Normal file
131
vtkc/Vertex.h
Normal file
@ -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 <vtkc/ErrorCode.h>
|
||||||
|
#include <vtkc/Shapes.h>
|
||||||
|
|
||||||
|
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<typename CoordType>
|
||||||
|
VTKC_EXEC inline vtkc::ErrorCode parametricCenter(Vertex, CoordType&&) noexcept
|
||||||
|
{
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
template<typename CoordType>
|
||||||
|
VTKC_EXEC inline vtkc::ErrorCode parametricPoint(Vertex, IdComponent pointId, CoordType&&) noexcept
|
||||||
|
{
|
||||||
|
if (pointId == 0)
|
||||||
|
{
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
return ErrorCode::INVALID_POINT_ID;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
template<typename CoordType>
|
||||||
|
VTKC_EXEC inline ComponentType<CoordType> parametricDistance(Vertex, const CoordType&) noexcept
|
||||||
|
{
|
||||||
|
return ComponentType<CoordType>{1};
|
||||||
|
}
|
||||||
|
|
||||||
|
template<typename CoordType>
|
||||||
|
VTKC_EXEC inline bool cellInside(Vertex, const CoordType&) noexcept
|
||||||
|
{
|
||||||
|
return false;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Values, typename CoordType, typename Result>
|
||||||
|
VTKC_EXEC inline vtkc::ErrorCode interpolate(
|
||||||
|
Vertex,
|
||||||
|
const Values& values,
|
||||||
|
const CoordType&,
|
||||||
|
Result&& result) noexcept
|
||||||
|
{
|
||||||
|
using T = ComponentType<Result>;
|
||||||
|
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
|
||||||
|
{
|
||||||
|
component(result, c) = static_cast<T>(values.getValue(0, c));
|
||||||
|
}
|
||||||
|
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Points, typename Values, typename CoordType, typename Result>
|
||||||
|
VTKC_EXEC inline vtkc::ErrorCode derivative(
|
||||||
|
Vertex,
|
||||||
|
const Points&,
|
||||||
|
const Values& values,
|
||||||
|
const CoordType&,
|
||||||
|
Result&& dx,
|
||||||
|
Result&& dy,
|
||||||
|
Result&& dz) noexcept
|
||||||
|
{
|
||||||
|
using T = ComponentType<Result>;
|
||||||
|
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
|
||||||
|
{
|
||||||
|
component(dx, c) = component(dy, c) = component(dz, c) = T{0};
|
||||||
|
}
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Points, typename PCoordType, typename WCoordType>
|
||||||
|
VTKC_EXEC inline vtkc::ErrorCode parametricToWorld(
|
||||||
|
Vertex,
|
||||||
|
const Points& points,
|
||||||
|
const PCoordType&,
|
||||||
|
WCoordType&& wcoords) noexcept
|
||||||
|
{
|
||||||
|
using T = ComponentType<WCoordType>;
|
||||||
|
component(wcoords, 0) = static_cast<T>(points.getValue(0, 0));
|
||||||
|
component(wcoords, 1) = static_cast<T>(points.getValue(0, 1));
|
||||||
|
component(wcoords, 2) = static_cast<T>(points.getValue(0, 2));
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Points, typename WCoordType, typename PCoordType>
|
||||||
|
VTKC_EXEC inline vtkc::ErrorCode worldToParametric(
|
||||||
|
Vertex,
|
||||||
|
const Points&,
|
||||||
|
const WCoordType&,
|
||||||
|
PCoordType&&) noexcept
|
||||||
|
{
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
} //namespace vtkc
|
||||||
|
|
||||||
|
#endif //vtk_c_Vertex_h
|
136
vtkc/Voxel.h
Normal file
136
vtkc/Voxel.h
Normal file
@ -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 <vtkc/ErrorCode.h>
|
||||||
|
#include <vtkc/Hexahedron.h>
|
||||||
|
|
||||||
|
#include <vtkc/internal/Common.h>
|
||||||
|
|
||||||
|
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 <typename Points, typename T>
|
||||||
|
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<T>(points.getValue(6, i) - points.getValue(0, i));
|
||||||
|
if (spacing[i] == T(0))
|
||||||
|
{
|
||||||
|
zeros |= 1 << i;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return zeros;
|
||||||
|
}
|
||||||
|
|
||||||
|
} // internal
|
||||||
|
|
||||||
|
template <typename Points, typename Values, typename CoordType, typename Result>
|
||||||
|
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<typename Values::ValueType>;
|
||||||
|
using ResultCompType = ComponentType<Result>;
|
||||||
|
|
||||||
|
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<ResultCompType>(dvdp[0] / spacing[0]);
|
||||||
|
component(dy, c) = static_cast<ResultCompType>(dvdp[1] / spacing[1]);
|
||||||
|
component(dz, c) = static_cast<ResultCompType>(dvdp[2] / spacing[2]);
|
||||||
|
}
|
||||||
|
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Points, typename PCoordType, typename WCoordType>
|
||||||
|
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<WCoordType>;
|
||||||
|
|
||||||
|
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<T>(component(pcoords, 0)));
|
||||||
|
component(wcoords, 1) = points.getValue(0, 1) +
|
||||||
|
(spacing[1] * static_cast<T>(component(pcoords, 1)));
|
||||||
|
component(wcoords, 2) = points.getValue(0, 2) +
|
||||||
|
(spacing[2] * static_cast<T>(component(pcoords, 2)));
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Points, typename WCoordType, typename PCoordType>
|
||||||
|
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<PCoordType>;
|
||||||
|
|
||||||
|
T spacing[3];
|
||||||
|
if (internal::getVoxelSpacing(points, spacing) != 0)
|
||||||
|
{
|
||||||
|
return ErrorCode::DEGENERATE_CELL_DETECTED;
|
||||||
|
}
|
||||||
|
|
||||||
|
component(pcoords, 0) = static_cast<T>(component(wcoords, 0) - points.getValue(0, 0)) /
|
||||||
|
spacing[0];
|
||||||
|
component(pcoords, 1) = static_cast<T>(component(wcoords, 1) - points.getValue(0, 1)) /
|
||||||
|
spacing[1];
|
||||||
|
component(pcoords, 2) = static_cast<T>(component(wcoords, 2) - points.getValue(0, 2)) /
|
||||||
|
spacing[2];
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
} //namespace vtkc
|
||||||
|
|
||||||
|
#endif //vtk_c_Voxel_h
|
228
vtkc/Wedge.h
Normal file
228
vtkc/Wedge.h
Normal file
@ -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 <vtkc/ErrorCode.h>
|
||||||
|
#include <vtkc/Shapes.h>
|
||||||
|
|
||||||
|
#include <vtkc/internal/Common.h>
|
||||||
|
|
||||||
|
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<typename CoordType>
|
||||||
|
VTKC_EXEC inline vtkc::ErrorCode parametricCenter(Wedge, CoordType&& pcoords) noexcept
|
||||||
|
{
|
||||||
|
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
|
||||||
|
|
||||||
|
using T = ComponentType<CoordType>;
|
||||||
|
component(pcoords, 0) = T(1)/T(3);
|
||||||
|
component(pcoords, 1) = T(1)/T(3);
|
||||||
|
component(pcoords, 2) = 0.5f;
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
template<typename CoordType>
|
||||||
|
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<typename CoordType>
|
||||||
|
VTKC_EXEC inline ComponentType<CoordType> parametricDistance(Wedge, const CoordType& pcoords) noexcept
|
||||||
|
{
|
||||||
|
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
|
||||||
|
return internal::findParametricDistance(pcoords, 3);
|
||||||
|
}
|
||||||
|
|
||||||
|
template<typename CoordType>
|
||||||
|
VTKC_EXEC inline bool cellInside(Wedge, const CoordType& pcoords) noexcept
|
||||||
|
{
|
||||||
|
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
|
||||||
|
|
||||||
|
using T = ComponentType<CoordType>;
|
||||||
|
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 <typename Values, typename CoordType, typename Result>
|
||||||
|
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<typename Values::ValueType>;
|
||||||
|
|
||||||
|
auto p0 = static_cast<T>(component(pcoords, 0));
|
||||||
|
auto p1 = static_cast<T>(component(pcoords, 1));
|
||||||
|
auto p2 = static_cast<T>(component(pcoords, 2));
|
||||||
|
auto sm = T{1} - p0 - p1;
|
||||||
|
|
||||||
|
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
|
||||||
|
{
|
||||||
|
T t1 = static_cast<T>(values.getValue(0, c)) * sm +
|
||||||
|
static_cast<T>(values.getValue(1, c)) * p0 +
|
||||||
|
static_cast<T>(values.getValue(2, c)) * p1;
|
||||||
|
T t2 = static_cast<T>(values.getValue(3, c)) * sm +
|
||||||
|
static_cast<T>(values.getValue(4, c)) * p0 +
|
||||||
|
static_cast<T>(values.getValue(5, c)) * p1;
|
||||||
|
component(result, c) = static_cast<ComponentType<Result>>(internal::lerp(t1, t2, p2));
|
||||||
|
}
|
||||||
|
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
namespace internal
|
||||||
|
{
|
||||||
|
|
||||||
|
template <typename Values, typename CoordType, typename Result>
|
||||||
|
VTKC_EXEC inline void parametricDerivative(
|
||||||
|
Wedge, const Values& values, IdComponent comp, const CoordType& pcoords, Result&& result) noexcept
|
||||||
|
{
|
||||||
|
using T = internal::ClosestFloatType<typename Values::ValueType>;
|
||||||
|
auto p0 = static_cast<T>(component(pcoords, 0));
|
||||||
|
auto p1 = static_cast<T>(component(pcoords, 1));
|
||||||
|
auto p2 = static_cast<T>(component(pcoords, 2));
|
||||||
|
auto rm = T{1} - p2;
|
||||||
|
auto sm = T{1} - p0 - p1;
|
||||||
|
|
||||||
|
T dr = (static_cast<T>(values.getValue(0, comp)) * -rm) +
|
||||||
|
(static_cast<T>(values.getValue(1, comp)) * rm) +
|
||||||
|
(static_cast<T>(values.getValue(3, comp)) * -p2) +
|
||||||
|
(static_cast<T>(values.getValue(4, comp)) * p2);
|
||||||
|
|
||||||
|
T ds = (static_cast<T>(values.getValue(0, comp)) * -rm) +
|
||||||
|
(static_cast<T>(values.getValue(2, comp)) * rm) +
|
||||||
|
(static_cast<T>(values.getValue(3, comp)) * -p2) +
|
||||||
|
(static_cast<T>(values.getValue(5, comp)) * p2);
|
||||||
|
|
||||||
|
T dt = (static_cast<T>(values.getValue(0, comp)) * -sm) +
|
||||||
|
(static_cast<T>(values.getValue(1, comp)) * -p0) +
|
||||||
|
(static_cast<T>(values.getValue(2, comp)) * -p1) +
|
||||||
|
(static_cast<T>(values.getValue(3, comp)) * sm) +
|
||||||
|
(static_cast<T>(values.getValue(4, comp)) * p0) +
|
||||||
|
(static_cast<T>(values.getValue(5, comp)) * p1);
|
||||||
|
|
||||||
|
component(result, 0) = static_cast<ComponentType<Result>>(dr);
|
||||||
|
component(result, 1) = static_cast<ComponentType<Result>>(ds);
|
||||||
|
component(result, 2) = static_cast<ComponentType<Result>>(dt);
|
||||||
|
}
|
||||||
|
|
||||||
|
} // internal
|
||||||
|
|
||||||
|
template <typename Points, typename Values, typename CoordType, typename Result>
|
||||||
|
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<Result>(dx),
|
||||||
|
std::forward<Result>(dy),
|
||||||
|
std::forward<Result>(dz));
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Points, typename PCoordType, typename WCoordType>
|
||||||
|
VTKC_EXEC inline vtkc::ErrorCode parametricToWorld(
|
||||||
|
Wedge,
|
||||||
|
const Points& points,
|
||||||
|
const PCoordType& pcoords,
|
||||||
|
WCoordType&& wcoords) noexcept
|
||||||
|
{
|
||||||
|
return interpolate(Wedge{}, points, pcoords, std::forward<WCoordType>(wcoords));
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Points, typename WCoordType, typename PCoordType>
|
||||||
|
VTKC_EXEC inline vtkc::ErrorCode worldToParametric(
|
||||||
|
Wedge,
|
||||||
|
const Points& points,
|
||||||
|
const WCoordType& wcoords,
|
||||||
|
PCoordType&& pcoords) noexcept
|
||||||
|
{
|
||||||
|
return internal::worldToParametric3D(Wedge{}, points, wcoords, std::forward<PCoordType>(pcoords));
|
||||||
|
}
|
||||||
|
|
||||||
|
} // vtkc
|
||||||
|
|
||||||
|
#endif // vtk_c_Wedge_h
|
310
vtkc/internal/Common.h
Normal file
310
vtkc/internal/Common.h
Normal file
@ -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 <vtkc/FieldAccessor.h>
|
||||||
|
#include <vtkc/internal/Math.h>
|
||||||
|
|
||||||
|
#include <cstdint>
|
||||||
|
#include <type_traits>
|
||||||
|
|
||||||
|
#define VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType) \
|
||||||
|
static_assert(std::is_floating_point<ComponentType<PCoordType>>::value, \
|
||||||
|
"parametric coordinates should be of floating point type");
|
||||||
|
|
||||||
|
namespace vtkc
|
||||||
|
{
|
||||||
|
namespace internal
|
||||||
|
{
|
||||||
|
|
||||||
|
///=========================================================================
|
||||||
|
template<typename CoordType>
|
||||||
|
VTKC_EXEC inline ComponentType<CoordType> findParametricDistance(
|
||||||
|
const CoordType& pvals, IdComponent numVals) noexcept
|
||||||
|
{
|
||||||
|
using T = ComponentType<CoordType>;
|
||||||
|
|
||||||
|
T pDistMax = 0.0f;
|
||||||
|
for (IdComponent i = 0; i < numVals; ++i)
|
||||||
|
{
|
||||||
|
ComponentType<CoordType> 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 <typename Values, typename CoordType, typename Result> \
|
||||||
|
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 <typename T>
|
||||||
|
class Space2D
|
||||||
|
{
|
||||||
|
public:
|
||||||
|
explicit VTKC_EXEC Space2D(const Vector<T, 3>& origin, const Vector<T, 3>& p1, const Vector<T, 3>& 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<T, 2> to2DPoint(Vector<T, 3> pt) const noexcept
|
||||||
|
{
|
||||||
|
pt -= this->Origin;
|
||||||
|
return Vector<T, 2>{ internal::dot(pt, this->XAxis), internal::dot(pt, this->YAxis) };
|
||||||
|
}
|
||||||
|
|
||||||
|
VTKC_EXEC Vector<T, 3> to3DVec(const Vector<T, 2>& vec) const noexcept
|
||||||
|
{
|
||||||
|
return (this->XAxis * vec[0]) + (this->YAxis * vec[1]);
|
||||||
|
}
|
||||||
|
|
||||||
|
private:
|
||||||
|
Vector<T, 3> Origin;
|
||||||
|
Vector<T, 3> XAxis, YAxis;
|
||||||
|
};
|
||||||
|
|
||||||
|
template <typename CellTag, typename Points, typename PCoords, typename T>
|
||||||
|
VTKC_EXEC inline void jacobian2D(
|
||||||
|
CellTag tag, const Points& points, const PCoords& pcoords, Matrix<T, 2, 2>& 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 <typename CellTag, typename Points, typename Values, typename CoordType, typename Result>
|
||||||
|
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<typename Values::ValueType>;
|
||||||
|
using ResultCompType = ComponentType<Result>;
|
||||||
|
|
||||||
|
constexpr IdComponent numPoints = CellTag{}.numberOfPoints();
|
||||||
|
|
||||||
|
Vector<ProcessingType, 3> pts[numPoints];
|
||||||
|
for (int i = 0; i < numPoints; ++i)
|
||||||
|
{
|
||||||
|
points.getTuple(i, pts[i]);
|
||||||
|
}
|
||||||
|
|
||||||
|
// 2-D coordinate system on the cell's plane
|
||||||
|
Space2D<ProcessingType> planeSpace(pts[0], pts[1], pts[numPoints - 1]);
|
||||||
|
Vector<ProcessingType, 2> pts2d[numPoints];
|
||||||
|
for (int i = 0; i < numPoints; ++i)
|
||||||
|
{
|
||||||
|
pts2d[i] = planeSpace.to2DPoint(pts[i]);
|
||||||
|
}
|
||||||
|
|
||||||
|
Matrix<ProcessingType, 2, 2> jacobian;
|
||||||
|
jacobian2D(tag, makeFieldAccessorNestedSOA(pts2d, 2), pcoords, jacobian);
|
||||||
|
Matrix<ProcessingType, 2, 2> invJacobian;
|
||||||
|
VTKC_RETURN_ON_ERROR(matrixInverse(jacobian, invJacobian))
|
||||||
|
|
||||||
|
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
|
||||||
|
{
|
||||||
|
Vector<ProcessingType, 2> dvdp;
|
||||||
|
parametricDerivative(tag, values, c, pcoords, dvdp);
|
||||||
|
auto d2D = matrixMultiply(dvdp, invJacobian);
|
||||||
|
auto d3D = planeSpace.to3DVec(d2D);
|
||||||
|
|
||||||
|
component(dx, c) = static_cast<ResultCompType>(d3D[0]);
|
||||||
|
component(dy, c) = static_cast<ResultCompType>(d3D[1]);
|
||||||
|
component(dz, c) = static_cast<ResultCompType>(d3D[2]);
|
||||||
|
}
|
||||||
|
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename CellTag, typename Points, typename WCoordType, typename PCoordType>
|
||||||
|
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<PCoordType>;
|
||||||
|
|
||||||
|
constexpr IdComponent numPoints = CellTag{}.numberOfPoints();
|
||||||
|
|
||||||
|
Vector<TIn, 3> pts[numPoints];
|
||||||
|
for (int i = 0; i < numPoints; ++i)
|
||||||
|
{
|
||||||
|
points.getTuple(i, pts[i]);
|
||||||
|
}
|
||||||
|
|
||||||
|
// 2-D coordinate system on the cell's plane
|
||||||
|
Space2D<TIn> planeSpace(pts[0], pts[1], pts[numPoints - 1]);
|
||||||
|
Vector<TIn, 2> pts2d[numPoints];
|
||||||
|
for (int i = 0; i < numPoints; ++i)
|
||||||
|
{
|
||||||
|
pts2d[i] = planeSpace.to2DPoint(pts[i]);
|
||||||
|
}
|
||||||
|
|
||||||
|
auto jacobianEvaluator = [&pts2d](const Vector<TOut, 2>& pc, Matrix<TIn, 2, 2>& jacobian) {
|
||||||
|
jacobian2D(CellTag{}, makeFieldAccessorNestedSOA(pts2d, 2), pc, jacobian);
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
};
|
||||||
|
|
||||||
|
auto functionEvaluator = [&points, &planeSpace](const Vector<TOut, 2>& pc, Vector<TIn, 2>& wc) {
|
||||||
|
Vector<TIn, 3> wc3(0);
|
||||||
|
VTKC_RETURN_ON_ERROR(parametricToWorld(CellTag{}, points, pc, wc3))
|
||||||
|
wc = planeSpace.to2DPoint(wc3);
|
||||||
|
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
};
|
||||||
|
|
||||||
|
Vector<TIn, 3> wcVec{component(wcoords, 0), component(wcoords, 1), component(wcoords, 2)};
|
||||||
|
Vector<TOut, 2> 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 <typename CellTag, typename Points, typename PCoords, typename T>
|
||||||
|
VTKC_EXEC inline void jacobian3D(
|
||||||
|
CellTag tag, const Points& points, const PCoords& pcoords, Matrix<T, 3, 3>& 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 <typename CellTag, typename Points, typename Values, typename CoordType, typename Result>
|
||||||
|
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<typename Values::ValueType>;
|
||||||
|
using ResultCompType = ComponentType<Result>;
|
||||||
|
|
||||||
|
Matrix<ProcessingType, 3, 3> jacobian;
|
||||||
|
jacobian3D(tag, points, pcoords, jacobian);
|
||||||
|
Matrix<ProcessingType, 3, 3> invJacobian;
|
||||||
|
VTKC_RETURN_ON_ERROR(matrixInverse(jacobian, invJacobian))
|
||||||
|
|
||||||
|
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
|
||||||
|
{
|
||||||
|
Vector<ProcessingType, 3> dvdp;
|
||||||
|
parametricDerivative(tag, values, c, pcoords, dvdp);
|
||||||
|
auto deriv = matrixMultiply(dvdp, invJacobian);
|
||||||
|
component(dx, c) = static_cast<ResultCompType>(deriv[0]);
|
||||||
|
component(dy, c) = static_cast<ResultCompType>(deriv[1]);
|
||||||
|
component(dz, c) = static_cast<ResultCompType>(deriv[2]);
|
||||||
|
}
|
||||||
|
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename CellTag, typename Points, typename WCoordType, typename PCoordType>
|
||||||
|
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<PCoordType>;
|
||||||
|
|
||||||
|
auto jacobianEvaluator = [tag, &points](const Vector<TOut, 3>& pc, Matrix<TIn, 3, 3>& jacobian) {
|
||||||
|
jacobian3D(tag, points, pc, jacobian);
|
||||||
|
return ErrorCode::SUCCESS;
|
||||||
|
};
|
||||||
|
|
||||||
|
auto functionEvaluator = [tag, &points](const Vector<TOut, 3>& pc, Vector<TIn, 3>& wc) {
|
||||||
|
return parametricToWorld(tag, points, pc, wc);
|
||||||
|
};
|
||||||
|
|
||||||
|
internal::Vector<TIn, 3> wcVec{component(wcoords, 0), component(wcoords, 1), component(wcoords, 2)};
|
||||||
|
internal::Vector<TOut, 3> 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
|
38
vtkc/internal/Config.h
Normal file
38
vtkc/internal/Config.h
Normal file
@ -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 <cstdint>
|
||||||
|
#include <type_traits>
|
||||||
|
|
||||||
|
#ifdef __CUDACC__
|
||||||
|
# define VTKC_EXEC __device__ __host__
|
||||||
|
#else
|
||||||
|
# define VTKC_EXEC
|
||||||
|
#endif
|
||||||
|
|
||||||
|
namespace vtkc
|
||||||
|
{
|
||||||
|
|
||||||
|
namespace internal
|
||||||
|
{
|
||||||
|
template <typename T>
|
||||||
|
using ClosestFloatType =
|
||||||
|
typename std::enable_if<std::is_arithmetic<T>::value,
|
||||||
|
typename std::conditional<sizeof(T) <= 4, float, double>::type>::type;
|
||||||
|
}
|
||||||
|
|
||||||
|
using IdShape = std::int8_t;
|
||||||
|
using IdComponent = std::int32_t;
|
||||||
|
|
||||||
|
} // namespace vtkc
|
||||||
|
|
||||||
|
#endif // vtk_c_internal_Config_h
|
712
vtkc/internal/Math.h
Normal file
712
vtkc/internal/Math.h
Normal file
@ -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 <cmath>
|
||||||
|
|
||||||
|
#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 <typename T, IdComponent Dim>
|
||||||
|
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 <typename... Ts>
|
||||||
|
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 <typename T, IdComponent Dim>
|
||||||
|
VTKC_EXEC
|
||||||
|
Vector<T, Dim>& operator+=(Vector<T, Dim>& v1, const Vector<T, Dim>& v2) noexcept
|
||||||
|
{
|
||||||
|
for (IdComponent i = 0; i < Dim; ++i)
|
||||||
|
{
|
||||||
|
v1[i] += v2[i];
|
||||||
|
}
|
||||||
|
return v1;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T, IdComponent Dim>
|
||||||
|
VTKC_EXEC
|
||||||
|
Vector<T, Dim>& operator-=(Vector<T, Dim>& v1, const Vector<T, Dim>& v2) noexcept
|
||||||
|
{
|
||||||
|
for (IdComponent i = 0; i < Dim; ++i)
|
||||||
|
{
|
||||||
|
v1[i] -= v2[i];
|
||||||
|
}
|
||||||
|
return v1;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T, IdComponent Dim>
|
||||||
|
VTKC_EXEC
|
||||||
|
Vector<T, Dim>& operator*=(Vector<T, Dim>& v1, const Vector<T, Dim>& v2) noexcept
|
||||||
|
{
|
||||||
|
for (IdComponent i = 0; i < Dim; ++i)
|
||||||
|
{
|
||||||
|
v1[i] *= v2[i];
|
||||||
|
}
|
||||||
|
return v1;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T, IdComponent Dim>
|
||||||
|
VTKC_EXEC
|
||||||
|
Vector<T, Dim>& operator*=(Vector<T, Dim>& v, const T& s) noexcept
|
||||||
|
{
|
||||||
|
for (IdComponent i = 0; i < Dim; ++i)
|
||||||
|
{
|
||||||
|
v[i] *= s;
|
||||||
|
}
|
||||||
|
return v;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T, IdComponent Dim>
|
||||||
|
VTKC_EXEC
|
||||||
|
Vector<T, Dim>& operator/=(Vector<T, Dim>& v1, const Vector<T, Dim>& v2) noexcept
|
||||||
|
{
|
||||||
|
for (IdComponent i = 0; i < Dim; ++i)
|
||||||
|
{
|
||||||
|
v1[i] /= v2[i];
|
||||||
|
}
|
||||||
|
return v1;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T, IdComponent Dim>
|
||||||
|
VTKC_EXEC
|
||||||
|
Vector<T, Dim>& operator/=(Vector<T, Dim>& v, const T& s) noexcept
|
||||||
|
{
|
||||||
|
for (IdComponent i = 0; i < Dim; ++i)
|
||||||
|
{
|
||||||
|
v[i] /= s;
|
||||||
|
}
|
||||||
|
return v;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T, IdComponent Dim>
|
||||||
|
VTKC_EXEC
|
||||||
|
Vector<T, Dim> operator+(Vector<T, Dim> v1, const Vector<T, Dim>& v2) noexcept
|
||||||
|
{
|
||||||
|
return v1 += v2;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T, IdComponent Dim>
|
||||||
|
VTKC_EXEC
|
||||||
|
Vector<T, Dim> operator-(Vector<T, Dim> v1, const Vector<T, Dim>& v2) noexcept
|
||||||
|
{
|
||||||
|
return v1 -= v2;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T, IdComponent Dim>
|
||||||
|
VTKC_EXEC
|
||||||
|
Vector<T, Dim> operator*(Vector<T, Dim> v1, const Vector<T, Dim>& v2) noexcept
|
||||||
|
{
|
||||||
|
return v1 *= v2;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T, IdComponent Dim>
|
||||||
|
VTKC_EXEC
|
||||||
|
Vector<T, Dim> operator*(Vector<T, Dim> v, const T& s) noexcept
|
||||||
|
{
|
||||||
|
return v *= s;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T, IdComponent Dim>
|
||||||
|
VTKC_EXEC
|
||||||
|
Vector<T, Dim> operator*(const T& s, Vector<T, Dim> v) noexcept
|
||||||
|
{
|
||||||
|
return v *= s;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T, IdComponent Dim>
|
||||||
|
VTKC_EXEC
|
||||||
|
Vector<T, Dim> operator/(Vector<T, Dim> v1, const Vector<T, Dim>& v2) noexcept
|
||||||
|
{
|
||||||
|
return v1 /= v2;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T, IdComponent Dim>
|
||||||
|
VTKC_EXEC
|
||||||
|
Vector<T, Dim> operator/(Vector<T, Dim> v, const T& s) noexcept
|
||||||
|
{
|
||||||
|
return v /= s;
|
||||||
|
}
|
||||||
|
|
||||||
|
//---------------------------------------------------------------------------
|
||||||
|
template <typename T, IdComponent Dim>
|
||||||
|
VTKC_EXEC
|
||||||
|
T dot(const Vector<T, Dim>& v1, const Vector<T, Dim>& v2) noexcept
|
||||||
|
{
|
||||||
|
T result{};
|
||||||
|
for (IdComponent i = 0; i < Dim; ++i)
|
||||||
|
{
|
||||||
|
result += v1[i] * v2[i];
|
||||||
|
}
|
||||||
|
return result;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T>
|
||||||
|
VTKC_EXEC
|
||||||
|
Vector<T, 3> cross(const Vector<T, 3>& v1, const Vector<T, 3>& v2) noexcept
|
||||||
|
{
|
||||||
|
return Vector<T, 3>(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 <typename T, IdComponent Dim>
|
||||||
|
VTKC_EXEC
|
||||||
|
T magnitude(const Vector<T, Dim>& v)
|
||||||
|
{
|
||||||
|
return VTKC_MATH_CALL(sqrt, (dot(v, v)));
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T, IdComponent Dim>
|
||||||
|
VTKC_EXEC
|
||||||
|
void normalize(Vector<T, Dim>& v)
|
||||||
|
{
|
||||||
|
v /= magnitude(v);
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T, IdComponent Dim>
|
||||||
|
VTKC_EXEC
|
||||||
|
Vector<T, Dim> normal(Vector<T, Dim> v)
|
||||||
|
{
|
||||||
|
normalize(v);
|
||||||
|
return v;
|
||||||
|
}
|
||||||
|
|
||||||
|
///=========================================================================
|
||||||
|
/// Basic matrix class and functions for internal use
|
||||||
|
|
||||||
|
template <typename T, IdComponent NumberOfRows, IdComponent NumberOfColumns>
|
||||||
|
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<T, NumberOfRows>& getColumn(IdComponent c) const
|
||||||
|
{
|
||||||
|
return this->Columns[c];
|
||||||
|
}
|
||||||
|
|
||||||
|
VTKC_EXEC
|
||||||
|
Vector<T, NumberOfColumns> getRow(IdComponent r) const
|
||||||
|
{
|
||||||
|
Vector<T, NumberOfColumns> 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<T, NumberOfRows> Columns[NumberOfColumns];
|
||||||
|
};
|
||||||
|
|
||||||
|
//---------------------------------------------------------------------------
|
||||||
|
template <typename T, IdComponent NRows, IdComponent NCols>
|
||||||
|
VTKC_EXEC
|
||||||
|
Matrix<T, NRows, NCols> operator+(Matrix<T, NRows, NCols> m1,
|
||||||
|
const Matrix<T, NRows, NCols>& m2) noexcept
|
||||||
|
{
|
||||||
|
return m1 += m2;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T, IdComponent NRows, IdComponent NCols>
|
||||||
|
VTKC_EXEC
|
||||||
|
Matrix<T, NRows, NCols> operator-(Matrix<T, NRows, NCols> m1,
|
||||||
|
const Matrix<T, NRows, NCols>& m2) noexcept
|
||||||
|
{
|
||||||
|
return m1 -= m2;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T, IdComponent NRows, IdComponent NCols>
|
||||||
|
VTKC_EXEC
|
||||||
|
Matrix<T, NRows, NCols> operator*(Matrix<T, NRows, NCols> m, const T& s) noexcept
|
||||||
|
{
|
||||||
|
return m *= s;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T, IdComponent NRows, IdComponent NCols>
|
||||||
|
VTKC_EXEC
|
||||||
|
Matrix<T, NRows, NCols> operator/(Matrix<T, NRows, NCols> m, const T& s) noexcept
|
||||||
|
{
|
||||||
|
return m /= s;
|
||||||
|
}
|
||||||
|
|
||||||
|
//---------------------------------------------------------------------------
|
||||||
|
template <typename T,
|
||||||
|
IdComponent NumRow,
|
||||||
|
IdComponent NumCol,
|
||||||
|
IdComponent NumInternal>
|
||||||
|
VTKC_EXEC
|
||||||
|
Matrix<T, NumRow, NumCol> matrixMultiply(
|
||||||
|
const Matrix<T, NumRow, NumInternal>& leftFactor,
|
||||||
|
const Matrix<T, NumInternal, NumCol>& rightFactor) noexcept
|
||||||
|
{
|
||||||
|
Matrix<T, NumRow, NumCol> 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 <typename T, IdComponent NumRow, IdComponent NumCol>
|
||||||
|
VTKC_EXEC
|
||||||
|
Vector<T, NumRow> matrixMultiply(
|
||||||
|
const Matrix<T, NumRow, NumCol>& leftFactor,
|
||||||
|
const Vector<T, NumCol>& rightFactor) noexcept
|
||||||
|
{
|
||||||
|
Vector<T, NumRow> result;
|
||||||
|
for (IdComponent rowIndex = 0; rowIndex < NumRow; ++rowIndex)
|
||||||
|
{
|
||||||
|
result[rowIndex] = dot(leftFactor.getRow(rowIndex), rightFactor);
|
||||||
|
}
|
||||||
|
return result;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T, IdComponent NumRow, IdComponent NumCol>
|
||||||
|
VTKC_EXEC
|
||||||
|
Vector<T, NumCol> matrixMultiply(
|
||||||
|
const Vector<T, NumRow>& leftFactor,
|
||||||
|
const Matrix<T, NumRow, NumCol>& rightFactor) noexcept
|
||||||
|
{
|
||||||
|
Vector<T, NumCol> 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 <typename T, IdComponent Size>
|
||||||
|
VTKC_EXEC
|
||||||
|
vtkc::ErrorCode matrixLUPFactorFindPivot(Matrix<T, Size, Size>& m,
|
||||||
|
Vector<IdComponent, Size>& permutation,
|
||||||
|
IdComponent topCornerIndex,
|
||||||
|
T& inversionParity) noexcept
|
||||||
|
{
|
||||||
|
constexpr T epsilon = std::is_floating_point<T>::value ?
|
||||||
|
(std::is_same<T, float>::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 <typename T, IdComponent Size>
|
||||||
|
VTKC_EXEC
|
||||||
|
void matrixLUPFactorFindUpperTriangleElements(
|
||||||
|
Matrix<T, Size, Size>& 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 <typename T, IdComponent Size>
|
||||||
|
VTKC_EXEC
|
||||||
|
vtkc::ErrorCode matrixLUPFactor(Matrix<T, Size, Size>& m,
|
||||||
|
Vector<IdComponent, Size>& 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 <typename T, IdComponent Size>
|
||||||
|
VTKC_EXEC
|
||||||
|
Vector<T, Size> matrixLUPSolve(
|
||||||
|
const Matrix<T, Size, Size>& LU,
|
||||||
|
const Vector<IdComponent, Size>& permutation,
|
||||||
|
const Vector<T, Size>& 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<T, Size> 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<T, Size> 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 <typename T, IdComponent Size>
|
||||||
|
VTKC_EXEC
|
||||||
|
vtkc::ErrorCode solveLinearSystem(const Matrix<T, Size, Size>& A,
|
||||||
|
const Vector<T, Size>& b,
|
||||||
|
Vector<T, Size>& x) noexcept
|
||||||
|
{
|
||||||
|
// First, we will make an LUP-factorization to help us.
|
||||||
|
Matrix<T, Size, Size> LU = A;
|
||||||
|
Vector<IdComponent, Size> 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 <typename T, IdComponent Size>
|
||||||
|
VTKC_EXEC
|
||||||
|
vtkc::ErrorCode matrixInverse(const Matrix<T, Size, Size>& A, Matrix<T, Size, Size>& invA) noexcept
|
||||||
|
{
|
||||||
|
// First, we will make an LUP-factorization to help us.
|
||||||
|
Matrix<T, Size, Size> LU = A;
|
||||||
|
Vector<IdComponent, Size> 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<T, Size> ICol(T(0));
|
||||||
|
for (IdComponent colIndex = 0; colIndex < Size; colIndex++)
|
||||||
|
{
|
||||||
|
ICol[colIndex] = 1;
|
||||||
|
Vector<T, Size> 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 <typename JacobianFunctor, typename FunctionFunctor, typename T, IdComponent Size>
|
||||||
|
VTKC_EXEC
|
||||||
|
inline vtkc::ErrorCode newtonsMethod(
|
||||||
|
const JacobianFunctor& jacobianEvaluator,
|
||||||
|
const FunctionFunctor& functionEvaluator,
|
||||||
|
const Vector<T, Size>& rhs,
|
||||||
|
Vector<T, Size>& 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<T, Size, Size> jacobian;
|
||||||
|
Vector<T, Size> fx(0);
|
||||||
|
|
||||||
|
VTKC_RETURN_ON_ERROR(jacobianEvaluator(x, jacobian))
|
||||||
|
VTKC_RETURN_ON_ERROR(functionEvaluator(x, fx))
|
||||||
|
|
||||||
|
Vector<T, Size> 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 <typename T>
|
||||||
|
VTKC_EXEC inline T lerp(T v0, T v1, T t)
|
||||||
|
{
|
||||||
|
static_assert(std::is_floating_point<T>::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
|
251
vtkc/vtkc.h
Normal file
251
vtkc/vtkc.h
Normal file
@ -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 <vtkc/Hexahedron.h>
|
||||||
|
#include <vtkc/Line.h>
|
||||||
|
#include <vtkc/Pixel.h>
|
||||||
|
#include <vtkc/Polygon.h>
|
||||||
|
#include <vtkc/Pyramid.h>
|
||||||
|
#include <vtkc/Quad.h>
|
||||||
|
#include <vtkc/Tetra.h>
|
||||||
|
#include <vtkc/Triangle.h>
|
||||||
|
#include <vtkc/Vertex.h>
|
||||||
|
#include <vtkc/Voxel.h>
|
||||||
|
#include <vtkc/Wedge.h>
|
||||||
|
|
||||||
|
#include <utility>
|
||||||
|
|
||||||
|
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<typename CoordType>
|
||||||
|
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<CoordType>(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<typename CoordType>
|
||||||
|
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<CoordType>(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<typename CoordType>
|
||||||
|
VTKC_EXEC inline ComponentType<CoordType> parametricDistance(Cell tag, const CoordType& pcoords) noexcept
|
||||||
|
{
|
||||||
|
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
|
||||||
|
|
||||||
|
ComponentType<CoordType> 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<typename CoordType>
|
||||||
|
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 <typename Values, typename CoordType, typename Result>
|
||||||
|
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>(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 <typename Points, typename Values, typename CoordType, typename Result>
|
||||||
|
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<Result>(dx),
|
||||||
|
std::forward<Result>(dy),
|
||||||
|
std::forward<Result>(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 <typename Points, typename PCoordType, typename WCoordType>
|
||||||
|
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<WCoordType>(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 <typename Points, typename WCoordType, typename PCoordType>
|
||||||
|
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<PCoordType>(pcoords)));
|
||||||
|
default:
|
||||||
|
status = ErrorCode::INVALID_SHAPE_ID;
|
||||||
|
}
|
||||||
|
return status;
|
||||||
|
}
|
||||||
|
|
||||||
|
} //namespace vtkc
|
||||||
|
|
||||||
|
#endif //vtk_c_vtkc_h
|
Loading…
Reference in New Issue
Block a user