diff --git a/docs/changelog/lagrangian-coherent-structures.md b/docs/changelog/lagrangian-coherent-structures.md new file mode 100644 index 000000000..8931252c9 --- /dev/null +++ b/docs/changelog/lagrangian-coherent-structures.md @@ -0,0 +1,13 @@ +# Lagrangian Coherent Strucutres (LCS) Filter for VTK-m + +The new filter `vtkm::filter::LagrangianStructures` is meant for Finite Time +Lyapunov Exponent (FTLE) calculation using VTK-m. +The filter allows users to calculate FTLE in two ways +1. Provide a dataset with a vector field, which will be used to generate a flow + map. +2. Provide a dataset containing a flow map, which can be readily used for the + FTLE field calculation. + +The filter returns a dataset with a point field named FTLE. +Is the input is strucutred and an auxiliary grid was not used, the filter will +add the field to the original dataset set, else a new structured dataset is returned. diff --git a/examples/lagrangian_structures/CMakeLists.txt b/examples/lagrangian_structures/CMakeLists.txt new file mode 100644 index 000000000..93d8351ea --- /dev/null +++ b/examples/lagrangian_structures/CMakeLists.txt @@ -0,0 +1,25 @@ +##============================================================================= +## +## Copyright (c) Kitware, Inc. +## All rights reserved. +## See LICENSE.txt for details. +## +## This software is distributed WITHOUT ANY WARRANTY; without even +## the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +## PURPOSE. See the above copyright notice for more information. +## +##============================================================================= +cmake_minimum_required(VERSION 3.8...3.14 FATAL_ERROR) +project(ParticleAdvection CXX) + +#Find the VTK-m package +find_package(VTKm REQUIRED QUIET) + +add_executable(ftle LagrangianStructures.cxx) +target_link_libraries(ftle PRIVATE vtkm_cont vtkm_worklet) +vtkm_add_target_information(ftle + MODIFY_CUDA_FLAGS + DEVICE_SOURCES LagrangianStructures.cxx) +if(TARGET vtkm::tbb) + target_compile_definitions(ftle PRIVATE BUILDING_TBB_VERSION) +endif() diff --git a/examples/lagrangian_structures/LagrangianStructures.cxx b/examples/lagrangian_structures/LagrangianStructures.cxx new file mode 100644 index 000000000..defaf8836 --- /dev/null +++ b/examples/lagrangian_structures/LagrangianStructures.cxx @@ -0,0 +1,58 @@ +//============================================================================= +// +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +// +//============================================================================= + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include + +int main(int argc, char** argv) +{ + vtkm::cont::Initialize(argc, argv); + + if (argc < 3) + { + std::cout << "Usage : flte " << std::endl; + } + std::string datasetName(argv[1]); + std::string variableName(argv[2]); + + std::cout << "Reading input dataset" << std::endl; + vtkm::cont::DataSet input; + vtkm::io::reader::VTKDataSetReader reader(datasetName); + input = reader.ReadDataSet(); + std::cout << "Read input dataset" << std::endl; + + vtkm::filter::LagrangianStructures lcsFilter; + lcsFilter.SetStepSize(0.025f); + lcsFilter.SetNumberOfSteps(500); + lcsFilter.SetAdvectionTime(0.025f * 500); + lcsFilter.SetOutputFieldName("gradient"); + lcsFilter.SetActiveField(variableName); + + vtkm::cont::DataSet output = lcsFilter.Execute(input); + vtkm::io::writer::VTKDataSetWriter writer("out.vtk"); + writer.WriteDataSet(output); + std::cout << "Written output dataset" << std::endl; + + return 0; +} diff --git a/vtkm/filter/CMakeLists.txt b/vtkm/filter/CMakeLists.txt index 265d0be01..cc5e090f8 100644 --- a/vtkm/filter/CMakeLists.txt +++ b/vtkm/filter/CMakeLists.txt @@ -41,6 +41,7 @@ set(headers Histogram.h ImageConnectivity.h Lagrangian.h + LagrangianStructures.h MarchingCubes.h Mask.h MaskPoints.h @@ -104,6 +105,7 @@ set(header_template_sources Histogram.hxx ImageConnectivity.hxx Lagrangian.hxx + LagrangianStructures.hxx MarchingCubes.hxx Mask.hxx MaskPoints.hxx diff --git a/vtkm/filter/LagrangianStructures.h b/vtkm/filter/LagrangianStructures.h new file mode 100644 index 000000000..c5df74b02 --- /dev/null +++ b/vtkm/filter/LagrangianStructures.h @@ -0,0 +1,92 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +//============================================================================ + +#ifndef vtk_m_filter_LagrangianStructures_h +#define vtk_m_filter_LagrangianStructures_h + +#include +#include +#include +#include + +namespace vtkm +{ +namespace filter +{ + +class LagrangianStructures : public vtkm::filter::FilterDataSetWithField +{ +public: + using SupportedTypes = vtkm::TypeListTagFieldVec3; + + using Scalar = vtkm::worklet::particleadvection::ScalarType; + using Vector = vtkm::Vec; + + LagrangianStructures(); + + void SetStepSize(Scalar s) { this->StepSize = s; } + Scalar GetStepSize() { return this->StepSize; } + + void SetNumberOfSteps(vtkm::Id n) { this->NumberOfSteps = n; } + vtkm::Id GetNumberOfSteps() { return this->NumberOfSteps; } + + void SetAdvectionTime(Scalar advectionTime) { this->AdvectionTime = advectionTime; } + Scalar GetAdvectionTime() { return this->AdvectionTime; } + + void SetUseAuxiliaryGrid(bool useAuxiliaryGrid) { this->UseAuxiliaryGrid = useAuxiliaryGrid; } + bool GetUseAuxiliaryGrid() { return this->UseAuxiliaryGrid; } + + void SetAuxiliaryGridDimensions(vtkm::Id3 auxiliaryDims) { this->AuxiliaryDims = auxiliaryDims; } + vtkm::Id3 GetAuxiliaryGridDimensions() { return this->AuxiliaryDims; } + + void SetUseFlowMapOutput(bool useFlowMapOutput) { this->UseFlowMapOutput = useFlowMapOutput; } + bool GetUseFlowMapOutput() { return this->UseFlowMapOutput; } + + void SetOutputFieldName(std::string outputFieldName) { this->OutputFieldName = outputFieldName; } + std::string GetOutputFieldName() { return this->OutputFieldName; } + + inline void SetFlowMapOutput(vtkm::cont::ArrayHandle& flowMap) + { + this->FlowMapOutput = flowMap; + } + inline vtkm::cont::ArrayHandle GetFlowMapOutput() { return this->FlowMapOutput; } + + template + VTKM_CONT vtkm::cont::DataSet DoExecute( + const vtkm::cont::DataSet& input, + const vtkm::cont::ArrayHandle, StorageType>& field, + const vtkm::filter::FieldMetadata& fieldMeta, + const vtkm::filter::PolicyBase& policy); + + //Map a new field onto the resulting dataset after running the filter + //this call is only valid + template + VTKM_CONT bool DoMapField(vtkm::cont::DataSet& result, + const vtkm::cont::ArrayHandle& input, + const vtkm::filter::FieldMetadata& fieldMeta, + vtkm::filter::PolicyBase policy); + +private: + Scalar StepSize; + vtkm::Id NumberOfSteps; + Scalar AdvectionTime; + bool UseAuxiliaryGrid = false; + vtkm::Id3 AuxiliaryDims; + bool UseFlowMapOutput = false; + std::string OutputFieldName; + vtkm::cont::ArrayHandle FlowMapOutput; +}; + +} // namespace filter +} // namespace vtkm + +#include + +#endif // vtk_m_filter_LagrangianStructures_h diff --git a/vtkm/filter/LagrangianStructures.hxx b/vtkm/filter/LagrangianStructures.hxx new file mode 100644 index 000000000..0ae7ad1f8 --- /dev/null +++ b/vtkm/filter/LagrangianStructures.hxx @@ -0,0 +1,146 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +//============================================================================ + +#include +#include +#include +#include +#include +#include + +#include + +namespace vtkm +{ +namespace filter +{ + +//----------------------------------------------------------------------------- +inline VTKM_CONT LagrangianStructures::LagrangianStructures() + : vtkm::filter::FilterDataSetWithField() +{ + OutputFieldName = std::string("FTLE"); +} + +//----------------------------------------------------------------------------- +template +inline VTKM_CONT vtkm::cont::DataSet LagrangianStructures::DoExecute( + const vtkm::cont::DataSet& input, + const vtkm::cont::ArrayHandle, StorageType>& field, + const vtkm::filter::FieldMetadata& fieldMeta, + const vtkm::filter::PolicyBase&) +{ + if (!fieldMeta.IsPointField()) + { + throw vtkm::cont::ErrorFilterExecution("Point field expected."); + } + + using Structured2DType = vtkm::cont::CellSetStructured<2>; + using Structured3DType = vtkm::cont::CellSetStructured<3>; + + using FieldHandle = vtkm::cont::ArrayHandle, StorageType>; + using GridEvaluator = vtkm::worklet::particleadvection::GridEvaluator; + using Integrator = vtkm::worklet::particleadvection::RK4Integrator; + + Scalar stepSize = this->GetStepSize(); + vtkm::Id numberOfSteps = this->GetNumberOfSteps(); + + vtkm::cont::CoordinateSystem coordinates = input.GetCoordinateSystem(); + vtkm::cont::DynamicCellSet cellset = input.GetCellSet(); + + vtkm::cont::DataSet lcsInput; + if (this->GetUseAuxiliaryGrid()) + { + vtkm::Id3 lcsGridDims = this->GetAuxiliaryGridDimensions(); + vtkm::Bounds inputBounds = coordinates.GetBounds(); + Vector origin(static_cast(inputBounds.X.Min), + static_cast(inputBounds.Y.Min), + static_cast(inputBounds.Z.Min)); + Vector spacing; + spacing[0] = + static_cast(inputBounds.X.Length()) / static_cast(lcsGridDims[0] - 1); + spacing[1] = + static_cast(inputBounds.Y.Length()) / static_cast(lcsGridDims[1] - 1); + spacing[2] = + static_cast(inputBounds.Z.Length()) / static_cast(lcsGridDims[2] - 1); + vtkm::cont::DataSetBuilderUniform uniformDatasetBuilder; + lcsInput = uniformDatasetBuilder.Create(lcsGridDims, origin, spacing); + } + else + { + // Check if input dataset is structured. + // If not, we cannot proceed. + if (!(cellset.IsType() || cellset.IsType())) + throw vtkm::cont::ErrorFilterExecution( + "Provided data is not structured, provide parameters for an auxiliary grid."); + lcsInput = input; + } + vtkm::cont::ArrayHandle lcsInputPoints, lcsOutputPoints; + vtkm::cont::ArrayCopy(lcsInput.GetCoordinateSystem().GetData(), lcsInputPoints); + if (this->GetUseFlowMapOutput()) + { + // Check if there is a 1:1 correspondense between the flow map + // and the input points + lcsOutputPoints = this->GetFlowMapOutput(); + if (lcsInputPoints.GetNumberOfValues() != lcsOutputPoints.GetNumberOfValues()) + throw vtkm::cont::ErrorFilterExecution( + "Provided flow map does not correspond to the input points for LCS filter."); + } + else + { + GridEvaluator evaluator(input.GetCoordinateSystem(), input.GetCellSet(), field); + Integrator integrator(evaluator, stepSize); + vtkm::worklet::ParticleAdvection particles; + vtkm::worklet::ParticleAdvectionResult advectionResult; + vtkm::cont::ArrayHandle advectionPoints; + vtkm::cont::ArrayCopy(lcsInputPoints, advectionPoints); + advectionResult = particles.Run(integrator, advectionPoints, numberOfSteps); + lcsOutputPoints = advectionResult.positions; + } + // FTLE output field + vtkm::cont::ArrayHandle outputField; + vtkm::FloatDefault advectionTime = this->GetAdvectionTime(); + + vtkm::cont::DynamicCellSet lcsCellSet = lcsInput.GetCellSet(); + if (lcsCellSet.IsType()) + { + using AnalysisType = vtkm::worklet::LagrangianStructures<2>; + AnalysisType ftleCalculator(advectionTime, lcsCellSet); + vtkm::worklet::DispatcherMapField dispatcher(ftleCalculator); + dispatcher.Invoke(lcsInputPoints, lcsOutputPoints, outputField); + } + else if (lcsCellSet.IsType()) + { + using AnalysisType = vtkm::worklet::LagrangianStructures<3>; + AnalysisType ftleCalculator(advectionTime, lcsCellSet); + vtkm::worklet::DispatcherMapField dispatcher(ftleCalculator); + dispatcher.Invoke(lcsInputPoints, lcsOutputPoints, outputField); + } + + vtkm::cont::DataSet output; + vtkm::cont::DataSetFieldAdd fieldAdder; + output.AddCoordinateSystem(lcsInput.GetCoordinateSystem()); + output.AddCellSet(lcsInput.GetCellSet()); + fieldAdder.AddPointField(output, this->GetOutputFieldName(), outputField); + return output; +} + +//----------------------------------------------------------------------------- +template +inline VTKM_CONT bool LagrangianStructures::DoMapField( + vtkm::cont::DataSet&, + const vtkm::cont::ArrayHandle&, + const vtkm::filter::FieldMetadata&, + vtkm::filter::PolicyBase) +{ + return false; +} +} +} // namespace vtkm::filter diff --git a/vtkm/filter/testing/CMakeLists.txt b/vtkm/filter/testing/CMakeLists.txt index a0823d1cf..d96ee4b51 100644 --- a/vtkm/filter/testing/CMakeLists.txt +++ b/vtkm/filter/testing/CMakeLists.txt @@ -34,6 +34,7 @@ set(unit_tests UnitTestHistogramFilter.cxx UnitTestImageConnectivityFilter.cxx UnitTestLagrangianFilter.cxx + UnitTestLagrangianStructuresFilter.cxx UnitTestMarchingCubesFilter.cxx UnitTestMaskFilter.cxx UnitTestMaskPointsFilter.cxx diff --git a/vtkm/filter/testing/UnitTestLagrangianStructuresFilter.cxx b/vtkm/filter/testing/UnitTestLagrangianStructuresFilter.cxx new file mode 100644 index 000000000..097a0204a --- /dev/null +++ b/vtkm/filter/testing/UnitTestLagrangianStructuresFilter.cxx @@ -0,0 +1,228 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +//============================================================================ + +#include +#include +#include +#include + +#include + +namespace auxiliary +{ +vtkm::FloatDefault vecData[200 * 3] = { + 0.000000f, 0.000000f, 0.000000f, -0.032369f, 0.000000f, 0.000000f, -0.061107f, 0.000000f, + 0.000000f, -0.083145f, 0.000000f, 0.000000f, -0.096136f, 0.000000f, 0.000000f, -0.098712f, + 0.000000f, 0.000000f, -0.090620f, 0.000000f, 0.000000f, -0.072751f, -0.000000f, 0.000000f, + -0.047041f, 0.000000f, 0.000000f, -0.016264f, 0.000000f, 0.000000f, 0.016263f, 0.000000f, + 0.000000f, 0.047040f, 0.000000f, 0.000000f, 0.072750f, 0.000000f, 0.000000f, 0.090620f, + 0.000000f, 0.000000f, 0.098712f, -0.000000f, 0.000000f, 0.096136f, 0.000000f, 0.000000f, + 0.083145f, -0.000000f, 0.000000f, 0.061108f, 0.000000f, 0.000000f, 0.032369f, 0.000000f, + 0.000000f, 0.000000f, 0.000000f, 0.000000f, 0.000000f, 0.106811f, 0.000000f, -0.030274f, + 0.100785f, 0.000000f, -0.057152f, 0.083925f, 0.000000f, -0.077763f, 0.058081f, 0.000000f, + -0.089914f, 0.026058f, 0.000000f, -0.092323f, -0.008686f, 0.000000f, -0.084755f, -0.042409f, + 0.000000f, -0.068042f, -0.071488f, 0.000000f, -0.043996f, -0.092800f, 0.000000f, -0.015211f, + -0.104060f, 0.000000f, 0.015210f, -0.104060f, 0.000000f, 0.043995f, -0.092801f, 0.000000f, + 0.068042f, -0.071489f, 0.000000f, 0.084754f, -0.042410f, 0.000000f, 0.092323f, -0.008687f, + 0.000000f, 0.089914f, 0.026057f, 0.000000f, 0.077763f, 0.058080f, 0.000000f, 0.057152f, + 0.083925f, 0.000000f, 0.030274f, 0.100785f, 0.000000f, 0.000000f, 0.106811f, 0.000000f, + 0.000000f, 0.200103f, 0.000000f, -0.024596f, 0.188813f, 0.000000f, -0.046433f, 0.157227f, + 0.000000f, -0.063178f, 0.108809f, 0.000000f, -0.073050f, 0.048817f, 0.000000f, -0.075007f, + -0.016272f, 0.000000f, -0.068858f, -0.079450f, 0.000000f, -0.055280f, -0.133927f, 0.000000f, + -0.035744f, -0.173854f, 0.000000f, -0.012358f, -0.194948f, 0.000000f, 0.012357f, -0.194949f, + 0.000000f, 0.035743f, -0.173855f, 0.000000f, 0.055280f, -0.133929f, 0.000000f, 0.068858f, + -0.079452f, 0.000000f, 0.075007f, -0.016274f, 0.000000f, 0.073050f, 0.048816f, 0.000000f, + 0.063178f, 0.108809f, 0.000000f, 0.046433f, 0.157227f, 0.000000f, 0.024596f, 0.188813f, + 0.000000f, 0.000000f, 0.200103f, 0.000000f, 0.000000f, 0.269034f, 0.000000f, -0.016018f, + 0.253856f, 0.000000f, -0.030240f, 0.211388f, 0.000000f, -0.041145f, 0.146292f, 0.000000f, + -0.047574f, 0.065633f, 0.000000f, -0.048849f, -0.021878f, 0.000000f, -0.044845f, -0.106820f, + 0.000000f, -0.036002f, -0.180062f, 0.000000f, -0.023279f, -0.233744f, 0.000000f, -0.008049f, + -0.262104f, 0.000000f, 0.008048f, -0.262105f, 0.000000f, 0.023278f, -0.233745f, 0.000000f, + 0.036001f, -0.180065f, 0.000000f, 0.044844f, -0.106822f, 0.000000f, 0.048849f, -0.021880f, + 0.000000f, 0.047574f, 0.065632f, 0.000000f, 0.041145f, 0.146291f, 0.000000f, 0.030240f, + 0.211388f, 0.000000f, 0.016018f, 0.253856f, 0.000000f, 0.000000f, 0.269034f, 0.000000f, + 0.000000f, 0.305617f, 0.000000f, -0.005557f, 0.288374f, 0.000000f, -0.010491f, 0.240132f, + 0.000000f, -0.014274f, 0.166185f, 0.000000f, -0.016504f, 0.074558f, 0.000000f, -0.016947f, + -0.024852f, 0.000000f, -0.015557f, -0.121345f, 0.000000f, -0.012490f, -0.204547f, 0.000000f, + -0.008076f, -0.265527f, 0.000000f, -0.002792f, -0.297744f, 0.000000f, 0.002792f, -0.297745f, + 0.000000f, 0.008076f, -0.265529f, 0.000000f, 0.012490f, -0.204549f, 0.000000f, 0.015557f, + -0.121348f, 0.000000f, 0.016947f, -0.024855f, 0.000000f, 0.016504f, 0.074557f, 0.000000f, + 0.014274f, 0.166184f, 0.000000f, 0.010491f, 0.240132f, 0.000000f, 0.005557f, 0.288374f, + 0.000000f, 0.000000f, 0.305617f, 0.000000f, 0.000000f, 0.305617f, 0.000000f, 0.005557f, + 0.288374f, 0.000000f, 0.010491f, 0.240132f, 0.000000f, 0.014274f, 0.166185f, 0.000000f, + 0.016504f, 0.074558f, 0.000000f, 0.016947f, -0.024852f, 0.000000f, 0.015557f, -0.121345f, + 0.000000f, 0.012490f, -0.204547f, 0.000000f, 0.008076f, -0.265527f, 0.000000f, 0.002792f, + -0.297744f, 0.000000f, -0.002792f, -0.297745f, 0.000000f, -0.008076f, -0.265529f, 0.000000f, + -0.012490f, -0.204549f, 0.000000f, -0.015557f, -0.121348f, 0.000000f, -0.016947f, -0.024855f, + 0.000000f, -0.016504f, 0.074557f, 0.000000f, -0.014274f, 0.166184f, 0.000000f, -0.010491f, + 0.240132f, 0.000000f, -0.005557f, 0.288374f, 0.000000f, -0.000000f, 0.305617f, 0.000000f, + 0.000000f, 0.269034f, 0.000000f, 0.016018f, 0.253856f, 0.000000f, 0.030240f, 0.211388f, + 0.000000f, 0.041145f, 0.146292f, 0.000000f, 0.047574f, 0.065633f, 0.000000f, 0.048849f, + -0.021878f, 0.000000f, 0.044845f, -0.106820f, 0.000000f, 0.036002f, -0.180062f, 0.000000f, + 0.023279f, -0.233744f, 0.000000f, 0.008049f, -0.262104f, 0.000000f, -0.008048f, -0.262105f, + 0.000000f, -0.023278f, -0.233745f, 0.000000f, -0.036001f, -0.180065f, 0.000000f, -0.044844f, + -0.106822f, 0.000000f, -0.048849f, -0.021880f, 0.000000f, -0.047574f, 0.065632f, 0.000000f, + -0.041145f, 0.146291f, 0.000000f, -0.030240f, 0.211388f, 0.000000f, -0.016018f, 0.253856f, + 0.000000f, -0.000000f, 0.269034f, 0.000000f, 0.000000f, 0.200103f, 0.000000f, 0.024596f, + 0.188813f, 0.000000f, 0.046433f, 0.157227f, 0.000000f, 0.063178f, 0.108809f, 0.000000f, + 0.073050f, 0.048817f, 0.000000f, 0.075007f, -0.016272f, 0.000000f, 0.068858f, -0.079450f, + 0.000000f, 0.055280f, -0.133927f, 0.000000f, 0.035744f, -0.173854f, 0.000000f, 0.012358f, + -0.194948f, 0.000000f, -0.012357f, -0.194949f, 0.000000f, -0.035743f, -0.173855f, 0.000000f, + -0.055280f, -0.133929f, 0.000000f, -0.068858f, -0.079452f, 0.000000f, -0.075007f, -0.016274f, + 0.000000f, -0.073050f, 0.048816f, 0.000000f, -0.063178f, 0.108809f, 0.000000f, -0.046433f, + 0.157227f, 0.000000f, -0.024596f, 0.188813f, 0.000000f, -0.000000f, 0.200103f, 0.000000f, + 0.000000f, 0.106811f, 0.000000f, 0.030274f, 0.100785f, 0.000000f, 0.057152f, 0.083925f, + 0.000000f, 0.077763f, 0.058081f, 0.000000f, 0.089914f, 0.026058f, 0.000000f, 0.092323f, + -0.008686f, 0.000000f, 0.084755f, -0.042409f, 0.000000f, 0.068042f, -0.071488f, 0.000000f, + 0.043996f, -0.092800f, 0.000000f, 0.015211f, -0.104060f, 0.000000f, -0.015210f, -0.104060f, + 0.000000f, -0.043995f, -0.092801f, 0.000000f, -0.068042f, -0.071489f, 0.000000f, -0.084754f, + -0.042410f, 0.000000f, -0.092323f, -0.008687f, 0.000000f, -0.089914f, 0.026057f, 0.000000f, + -0.077763f, 0.058080f, 0.000000f, -0.057152f, 0.083925f, 0.000000f, -0.030274f, 0.100785f, + 0.000000f, -0.000000f, 0.106811f, 0.000000f, 0.000000f, 0.000000f, 0.000000f, 0.032369f, + 0.000000f, 0.000000f, 0.061107f, 0.000000f, 0.000000f, 0.083145f, 0.000000f, 0.000000f, + 0.096136f, 0.000000f, 0.000000f, 0.098712f, -0.000000f, 0.000000f, 0.090620f, -0.000000f, + 0.000000f, 0.072751f, -0.000000f, 0.000000f, 0.047041f, -0.000000f, 0.000000f, 0.016264f, + -0.000000f, 0.000000f, -0.016263f, -0.000000f, 0.000000f, -0.047040f, -0.000000f, 0.000000f, + -0.072750f, -0.000000f, 0.000000f, -0.090620f, -0.000000f, 0.000000f, -0.098712f, -0.000000f, + 0.000000f, -0.096136f, 0.000000f, 0.000000f, -0.083145f, 0.000000f, 0.000000f, -0.061108f, + 0.000000f, 0.000000f, -0.032369f, 0.000000f, 0.000000f, -0.000000f, 0.000000f, 0.000000f +}; + +vtkm::FloatDefault visitData[200] = { + 0.175776f, 0.193068f, 0.184354f, 0.162709f, 0.156567f, 0.165079f, 0.183523f, 0.176101f, 0.182635f, + 0.150985f, 0.127186f, 0.173021f, 0.173493f, 0.160812f, 0.124465f, 0.110838f, 0.135879f, 0.184853f, + 0.193068f, 0.000007f, 0.157161f, 0.144306f, 0.108385f, 0.086895f, 0.073086f, 0.092054f, 0.109178f, + 0.115175f, 0.144648f, 0.179178f, 0.177300f, 0.143649f, 0.113350f, 0.097418f, 0.081223f, 0.011488f, + 0.067968f, 0.108383f, 0.161194f, 0.191896f, 0.165629f, 0.133326f, 0.127635f, 0.058606f, 0.025823f, + 0.059042f, 0.100954f, 0.089389f, 0.124448f, 0.174170f, 0.174172f, 0.124446f, 0.089387f, 0.100955f, + 0.059041f, 0.025821f, 0.058604f, 0.127634f, 0.106871f, 0.188216f, 0.169607f, 0.135035f, 0.109767f, + 0.052501f, 0.029647f, 0.016691f, 0.076807f, 0.086811f, 0.121697f, 0.172919f, 0.172921f, 0.121695f, + 0.086809f, 0.076809f, 0.016692f, 0.029643f, 0.052500f, 0.109766f, 0.100262f, 0.181406f, 0.172943f, + 0.133487f, 0.098476f, 0.063748f, 0.031703f, 0.007183f, 0.062060f, 0.078550f, 0.119979f, 0.172034f, + 0.172036f, 0.119977f, 0.078548f, 0.062062f, 0.007186f, 0.031700f, 0.063747f, 0.098476f, 0.099140f, + 0.175363f, 0.175794f, 0.132061f, 0.090747f, 0.073996f, 0.033685f, 0.007585f, 0.053561f, 0.068431f, + 0.120009f, 0.171078f, 0.171080f, 0.120008f, 0.068429f, 0.053560f, 0.007588f, 0.033684f, 0.073997f, + 0.090747f, 0.132058f, 0.176183f, 0.177344f, 0.132153f, 0.091405f, 0.085764f, 0.052348f, 0.010314f, + 0.042262f, 0.064794f, 0.124248f, 0.168788f, 0.168789f, 0.124248f, 0.064795f, 0.042261f, 0.010313f, + 0.052345f, 0.085766f, 0.091405f, 0.132151f, 0.177345f, 0.177002f, 0.134660f, 0.101316f, 0.096766f, + 0.089878f, 0.029814f, 0.039441f, 0.096341f, 0.128119f, 0.162818f, 0.162820f, 0.128117f, 0.096345f, + 0.039441f, 0.029812f, 0.089876f, 0.096768f, 0.101316f, 0.134658f, 0.177003f, 0.174898f, 0.151468f, + 0.135607f, 0.108288f, 0.090863f, 0.080745f, 0.075302f, 0.100106f, 0.124623f, 0.161293f, 0.161294f, + 0.124626f, 0.100106f, 0.075302f, 0.080745f, 0.090864f, 0.108288f, 0.135607f, 0.151465f, 0.174900f, + 0.000001f, 0.161986f, 0.181619f, 0.169865f, 0.145339f, 0.163338f, 0.157883f, 0.171822f, 0.192594f, + 0.186741f, 0.186741f, 0.192594f, 0.171821f, 0.157882f, 0.163337f, 0.145338f, 0.169865f, 0.181618f, + 0.161986f, 0.000001f +}; + +vtkm::FloatDefault diffData[200] = { + 0.000000f, 0.000000f, 0.000000f, 0.000000f, 0.000000f, 0.000000f, 0.002152f, 0.003656f, 0.021269f, + 0.020149f, 0.043948f, 0.011656f, 0.006264f, 0.020559f, 0.040614f, 0.045728f, 0.026829f, 0.000501f, + 0.000000f, 0.175769f, 0.000001f, 0.000000f, 0.000000f, 0.000000f, 0.000000f, 0.000000f, 0.000000f, + 0.014023f, 0.001336f, 0.003623f, 0.001744f, 0.002334f, 0.015848f, 0.011760f, 0.010831f, 0.061597f, + 0.018926f, 0.000001f, 0.016888f, 0.034735f, 0.000000f, 0.000000f, 0.000000f, 0.000000f, 0.000000f, + 0.000001f, 0.000000f, 0.000000f, 0.000000f, 0.000001f, 0.000002f, 0.000002f, 0.000001f, 0.000000f, + 0.000000f, 0.000001f, 0.000000f, 0.000000f, 0.026455f, 0.022587f, 0.000000f, 0.000000f, 0.000000f, + 0.000001f, 0.000000f, 0.000001f, 0.000001f, 0.000000f, 0.000001f, 0.000001f, 0.000002f, 0.000001f, + 0.000002f, 0.000002f, 0.000001f, 0.000004f, 0.000001f, 0.000001f, 0.034773f, 0.011799f, 0.000000f, + 0.000000f, 0.000000f, 0.000000f, 0.000001f, 0.000000f, 0.000000f, 0.000000f, 0.000000f, 0.000001f, + 0.000002f, 0.000002f, 0.000002f, 0.000000f, 0.000003f, 0.000001f, 0.000001f, 0.000002f, 0.034346f, + 0.002419f, 0.000000f, 0.000000f, 0.000000f, 0.000000f, 0.000001f, 0.000001f, 0.000000f, 0.000000f, + 0.000000f, 0.000000f, 0.000002f, 0.000001f, 0.000002f, 0.000002f, 0.000002f, 0.000002f, 0.000000f, + 0.000002f, 0.000002f, 0.000389f, 0.000000f, 0.000000f, 0.000000f, 0.000000f, 0.000000f, 0.000001f, + 0.000000f, 0.000000f, 0.000000f, 0.000001f, 0.000001f, 0.000001f, 0.000000f, 0.000003f, 0.000002f, + 0.000000f, 0.000001f, 0.000002f, 0.000001f, 0.000001f, 0.000000f, 0.000000f, 0.000000f, 0.000001f, + 0.000000f, 0.000001f, 0.000000f, 0.000001f, 0.000000f, 0.000000f, 0.000002f, 0.000002f, 0.000001f, + 0.000001f, 0.000000f, 0.000001f, 0.000001f, 0.000002f, 0.000001f, 0.000001f, 0.000000f, 0.011404f, + 0.009266f, 0.013854f, 0.008200f, 0.000001f, 0.000002f, 0.000005f, 0.000006f, 0.000000f, 0.000000f, + 0.000004f, 0.000006f, 0.000002f, 0.000001f, 0.008199f, 0.013853f, 0.009266f, 0.011401f, 0.000002f, + 0.173409f, 0.040531f, 0.007484f, 0.012952f, 0.028796f, 0.004284f, 0.000001f, 0.000003f, 0.000006f, + 0.000009f, 0.000008f, 0.000006f, 0.000004f, 0.000002f, 0.004284f, 0.028796f, 0.012952f, 0.007483f, + 0.040531f, 0.173409f +}; + +using Scalar = vtkm::FloatDefault; +using Field = vtkm::Vec; + +void ValidateLCSFilterResult(const vtkm::cont::ArrayHandle& vtkmOutput, + const vtkm::cont::ArrayHandle& visitOutput, + const vtkm::cont::ArrayHandle& expDiff) +{ + VTKM_TEST_ASSERT(vtkmOutput.GetNumberOfValues() == 200, "Wrong number of values"); + + auto vtkmPortal = vtkmOutput.GetPortalConstControl(); + auto visitPortal = visitOutput.GetPortalConstControl(); + auto diffPortal = expDiff.GetPortalConstControl(); + + for (vtkm::Id index = 0; index < 200; index++) + { + auto vtkm = vtkmPortal.Get(index); + auto visit = visitPortal.Get(index); + auto diff = vtkm::Abs(vtkm - visit); + auto exp = diffPortal.Get(index); + VTKM_TEST_ASSERT(diff <= (vtkm::Abs(exp) + vtkm::Epsilon()), + "Calculated o/p is not as expected"); + } +} +} + +void TestLagrangianStructures() +{ + using Scalar = vtkm::FloatDefault; + using Field = vtkm::Vec; + + /*Construct dataset and vector field*/ + vtkm::cont::DataSet inputData; + vtkm::Id2 dims(20, 10); + vtkm::Vec origin(0.0f, 0.0f); + vtkm::Vec spacing; + spacing[0] = 2.0f / static_cast(dims[0] - 1); + spacing[1] = 1.0f / static_cast(dims[1] - 1); + vtkm::cont::DataSetBuilderUniform dataBuilder; + inputData = dataBuilder.Create(dims, origin, spacing); + + std::vector diffVec; + std::vector visitVec; + std::vector fieldVec; + for (vtkm::Id index = 0; index < 200; index++) + { + vtkm::Id flatidx = index * 3; + Field field; + field[0] = auxiliary::vecData[flatidx]; + field[1] = auxiliary::vecData[++flatidx]; + field[2] = auxiliary::vecData[++flatidx]; + fieldVec.push_back(field); + Scalar diff = auxiliary::diffData[index]; + diffVec.push_back(diff); + Scalar visit = auxiliary::visitData[index]; + visitVec.push_back(visit); + } + vtkm::cont::ArrayHandle diffHandle = vtkm::cont::make_ArrayHandle(diffVec); + vtkm::cont::ArrayHandle visitHandle = vtkm::cont::make_ArrayHandle(visitVec); + vtkm::cont::ArrayHandle fieldHandle = vtkm::cont::make_ArrayHandle(fieldVec); + vtkm::cont::DataSetFieldAdd fieldAdder; + fieldAdder.AddPointField(inputData, "velocity", fieldHandle); + + vtkm::filter::LagrangianStructures lagrangianStructures; + lagrangianStructures.SetStepSize(0.025f); + lagrangianStructures.SetNumberOfSteps(500); + lagrangianStructures.SetAdvectionTime(0.025f * 500); + lagrangianStructures.SetAdvectionTime(0.025f * 500); + lagrangianStructures.SetActiveField("velocity"); + vtkm::cont::DataSet outputData = lagrangianStructures.Execute(inputData); + + vtkm::cont::ArrayHandle FTLEField; + outputData.GetField("FTLE").GetData().CopyTo(FTLEField); + auxiliary::ValidateLCSFilterResult(FTLEField, visitHandle, diffHandle); +} + +int UnitTestLagrangianStructuresFilter(int argc, char* argv[]) +{ + return vtkm::cont::testing::Testing::Run(TestLagrangianStructures, argc, argv); +} diff --git a/vtkm/worklet/CMakeLists.txt b/vtkm/worklet/CMakeLists.txt index 783c2d017..229fd4dac 100644 --- a/vtkm/worklet/CMakeLists.txt +++ b/vtkm/worklet/CMakeLists.txt @@ -36,6 +36,7 @@ set(headers KdTree3D.h KernelSplatter.h Keys.h + LagrangianStructures.h Magnitude.h MarchingCubes.h Mask.h @@ -125,6 +126,7 @@ add_subdirectory(contourtree_augmented) add_subdirectory(cosmotools) add_subdirectory(gradient) add_subdirectory(histogram) +add_subdirectory(lcs) add_subdirectory(splatkernels) add_subdirectory(spatialstructure) add_subdirectory(tetrahedralize) diff --git a/vtkm/worklet/LagrangianStructures.h b/vtkm/worklet/LagrangianStructures.h new file mode 100644 index 000000000..27a363104 --- /dev/null +++ b/vtkm/worklet/LagrangianStructures.h @@ -0,0 +1,196 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +//============================================================================ +#ifndef vtk_m_worklet_LagrangianStructures_h +#define vtk_m_worklet_LagrangianStructures_h + +#include +#include +#include +#include + +#include +#include + +namespace vtkm +{ +namespace worklet +{ + +template +class LagrangianStructures; + +template <> +class LagrangianStructures<2> : public vtkm::worklet::WorkletMapField +{ +public: + using Scalar = vtkm::FloatDefault; + + VTKM_CONT + LagrangianStructures(Scalar endTime, vtkm::cont::DynamicCellSet cellSet) + : EndTime(endTime) + , GridData(cellSet) + { + } + + using ControlSignature = void(WholeArrayIn, WholeArrayIn, FieldOut); + + using ExecutionSignature = void(WorkIndex, _1, _2, _3); + + template + VTKM_EXEC void operator()(const vtkm::Id index, + const PointArray& input, + const PointArray& output, + Scalar& outputField) const + { + using Point = typename PointArray::ValueType; + + const vtkm::Vec neighborIndices = this->GridData.GetNeighborIndices(index); + + // Calculate Stretching / Squeezing + Point xin1 = input.Get(neighborIndices[0]); + Point xin2 = input.Get(neighborIndices[1]); + Point yin1 = input.Get(neighborIndices[2]); + Point yin2 = input.Get(neighborIndices[3]); + + Scalar xDiff = 1.0f / (xin2[0] - xin1[0]); + Scalar yDiff = 1.0f / (yin2[1] - yin1[1]); + + Point xout1 = output.Get(neighborIndices[0]); + Point xout2 = output.Get(neighborIndices[1]); + Point yout1 = output.Get(neighborIndices[2]); + Point yout2 = output.Get(neighborIndices[3]); + + // Total X gradient w.r.t X, Y + Scalar f1x = (xout2[0] - xout1[0]) * xDiff; + Scalar f1y = (yout2[0] - yout1[0]) * yDiff; + + // Total Y gradient w.r.t X, Y + Scalar f2x = (xout2[1] - xout1[1]) * xDiff; + Scalar f2y = (yout2[1] - yout1[1]) * yDiff; + + vtkm::Matrix jacobian; + vtkm::MatrixSetRow(jacobian, 0, vtkm::Vec(f1x, f1y)); + vtkm::MatrixSetRow(jacobian, 1, vtkm::Vec(f2x, f2y)); + + detail::ComputeLeftCauchyGreenTensor(jacobian); + + vtkm::Vec eigenValues; + detail::Jacobi(jacobian, eigenValues); + + Scalar delta = eigenValues[0]; + // Check if we need to clamp these values + // Also provide options. + // 1. FTLE + // 2. FLLE + // 3. Eigen Values (Min/Max) + //Scalar delta = trace + sqrtr; + // Given endTime is in units where start time is 0, + // else do endTime-startTime + // return value for computation + outputField = vtkm::Log(delta) / (static_cast(2.0f) * EndTime); + } + +public: + // To calculate FTLE field + Scalar EndTime; + // To assist in calculation of indices + detail::GridMetaData GridData; +}; + +template <> +class LagrangianStructures<3> : public vtkm::worklet::WorkletMapField +{ +public: + using Scalar = vtkm::FloatDefault; + + VTKM_CONT + LagrangianStructures(Scalar endTime, vtkm::cont::DynamicCellSet cellSet) + : EndTime(endTime) + , GridData(cellSet) + { + } + + using ControlSignature = void(WholeArrayIn, WholeArrayIn, FieldOut); + + using ExecutionSignature = void(WorkIndex, _1, _2, _3); + + /* + * Point position arrays are the input and the output positions of the particle advection. + */ + template + VTKM_EXEC void operator()(const vtkm::Id index, + const PointArray& input, + const PointArray& output, + Scalar& outputField) const + { + using Point = typename PointArray::ValueType; + + const vtkm::Vec neighborIndices = this->GridData.GetNeighborIndices(index); + + Point xin1 = input.Get(neighborIndices[0]); + Point xin2 = input.Get(neighborIndices[1]); + Point yin1 = input.Get(neighborIndices[2]); + Point yin2 = input.Get(neighborIndices[3]); + Point zin1 = input.Get(neighborIndices[4]); + Point zin2 = input.Get(neighborIndices[5]); + + Scalar xDiff = 1.0f / (xin2[0] - xin1[0]); + Scalar yDiff = 1.0f / (yin2[1] - yin1[1]); + Scalar zDiff = 1.0f / (zin2[2] - zin1[2]); + + Point xout1 = output.Get(neighborIndices[0]); + Point xout2 = output.Get(neighborIndices[1]); + Point yout1 = output.Get(neighborIndices[2]); + Point yout2 = output.Get(neighborIndices[3]); + Point zout1 = output.Get(neighborIndices[4]); + Point zout2 = output.Get(neighborIndices[5]); + + // Total X gradient w.r.t X, Y, Z + Scalar f1x = (xout2[0] - xout1[0]) * xDiff; + Scalar f1y = (yout2[0] - yout1[0]) * yDiff; + Scalar f1z = (zout2[0] - zout1[0]) * zDiff; + + // Total Y gradient w.r.t X, Y, Z + Scalar f2x = (xout2[1] - xout1[1]) * xDiff; + Scalar f2y = (yout2[1] - yout1[1]) * yDiff; + Scalar f2z = (zout2[1] - zout1[1]) * zDiff; + + // Total Z gradient w.r.t X, Y, Z + Scalar f3x = (xout2[2] - xout1[2]) * xDiff; + Scalar f3y = (yout2[2] - yout1[2]) * yDiff; + Scalar f3z = (zout2[2] - zout1[2]) * zDiff; + + vtkm::Matrix jacobian; + vtkm::MatrixSetRow(jacobian, 0, vtkm::Vec(f1x, f1y, f1z)); + vtkm::MatrixSetRow(jacobian, 1, vtkm::Vec(f2x, f2y, f2z)); + vtkm::MatrixSetRow(jacobian, 2, vtkm::Vec(f3x, f3y, f3z)); + + detail::ComputeLeftCauchyGreenTensor(jacobian); + + vtkm::Vec eigenValues; + detail::Jacobi(jacobian, eigenValues); + + Scalar delta = eigenValues[0]; + // Given endTime is in units where start time is 0. else do endTime-startTime + // return value for ftle computation + outputField = vtkm::Log(delta) / (static_cast(2.0f) * EndTime); + } + +public: + // To calculate FTLE field + Scalar EndTime; + // To assist in calculation of indices + detail::GridMetaData GridData; +}; + +} // worklet +} // vtkm + +#endif //vtk_m_worklet_LagrangianStructures_h diff --git a/vtkm/worklet/lcs/CMakeLists.txt b/vtkm/worklet/lcs/CMakeLists.txt new file mode 100644 index 000000000..2426d7a59 --- /dev/null +++ b/vtkm/worklet/lcs/CMakeLists.txt @@ -0,0 +1,16 @@ +##============================================================================ +## Copyright (c) Kitware, Inc. +## All rights reserved. +## See LICENSE.txt for details. +## +## This software is distributed WITHOUT ANY WARRANTY; without even +## the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +## PURPOSE. See the above copyright notice for more information. +##============================================================================ + +set(headers + GridMetaData.h + LagrangianStructureHelpers.h + ) + +vtkm_declare_headers(${headers}) diff --git a/vtkm/worklet/lcs/GridMetaData.h b/vtkm/worklet/lcs/GridMetaData.h new file mode 100644 index 000000000..68c2f6ff9 --- /dev/null +++ b/vtkm/worklet/lcs/GridMetaData.h @@ -0,0 +1,93 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +//============================================================================ +#ifndef vtk_m_worklet_lcs_GridMetaData_h +#define vtk_m_worklet_lcs_GridMetaData_h + +namespace vtkm +{ +namespace worklet +{ +namespace detail +{ + +class GridMetaData +{ +public: + using Structured2DType = vtkm::cont::CellSetStructured<2>; + using Structured3DType = vtkm::cont::CellSetStructured<3>; + + VTKM_CONT + GridMetaData(const vtkm::cont::DynamicCellSet cellSet) + { + if (cellSet.IsType()) + { + this->cellSet2D = true; + vtkm::Id2 dims = + cellSet.Cast().GetSchedulingRange(vtkm::TopologyElementTagPoint()); + this->Dims = vtkm::Id3(dims[0], dims[1], 1); + } + else + { + this->cellSet2D = false; + this->Dims = + cellSet.Cast().GetSchedulingRange(vtkm::TopologyElementTagPoint()); + } + this->PlaneSize = Dims[0] * Dims[1]; + this->RowSize = Dims[0]; + } + + VTKM_EXEC + bool IsCellSet2D() const { return this->cellSet2D; } + + VTKM_EXEC + void GetLogicalIndex(const vtkm::Id index, vtkm::Id3& logicalIndex) const + { + logicalIndex[0] = index % Dims[0]; + logicalIndex[1] = (index / Dims[0]) % Dims[1]; + if (this->cellSet2D) + logicalIndex[2] = 0; + else + logicalIndex[2] = index / (Dims[0] * Dims[1]); + } + + VTKM_EXEC + const vtkm::Vec GetNeighborIndices(const vtkm::Id index) const + { + vtkm::Vec indices; + vtkm::Id3 logicalIndex; + GetLogicalIndex(index, logicalIndex); + + // For differentials w.r.t delta in x + indices[0] = (logicalIndex[0] == 0) ? index : index - 1; + indices[1] = (logicalIndex[0] == Dims[0] - 1) ? index : index + 1; + // For differentials w.r.t delta in y + indices[2] = (logicalIndex[1] == 0) ? index : index - RowSize; + indices[3] = (logicalIndex[1] == Dims[1] - 1) ? index : index + RowSize; + if (!this->cellSet2D) + { + // For differentials w.r.t delta in z + indices[4] = (logicalIndex[2] == 0) ? index : index - PlaneSize; + indices[5] = (logicalIndex[2] == Dims[2] - 1) ? index : index + PlaneSize; + } + return indices; + } + +private: + bool cellSet2D = false; + vtkm::Id3 Dims; + vtkm::Id PlaneSize; + vtkm::Id RowSize; +}; + +} // namespace detail +} // namespace worklet +} // namespace vtkm + +#endif //vtk_m_worklet_lcs_GridMetaData_h diff --git a/vtkm/worklet/lcs/LagrangianStructureHelpers.h b/vtkm/worklet/lcs/LagrangianStructureHelpers.h new file mode 100644 index 000000000..54321b276 --- /dev/null +++ b/vtkm/worklet/lcs/LagrangianStructureHelpers.h @@ -0,0 +1,209 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +//============================================================================ +#ifndef vtk_m_worklet_lcs_LagrangianStructureHelpers_h +#define vtk_m_worklet_lcs_LagrangianStructureHelpers_h + +#include +#include +#include + +namespace vtkm +{ +namespace worklet +{ +namespace detail +{ + +template +VTKM_EXEC_CONT void ComputeLeftCauchyGreenTensor(vtkm::Matrix& jacobian) +{ + vtkm::Vec j1 = vtkm::MatrixGetRow(jacobian, 0); + vtkm::Vec j2 = vtkm::MatrixGetRow(jacobian, 1); + + // Left Cauchy Green Tensor is J*J^T + // j1[0] j1[1] | j1[0] j2[0] + // j2[0] j2[1] | j1[1] j2[1] + + T a = j1[0] * j1[0] + j1[1] * j1[1]; + T b = j1[0] * j2[0] + j1[1] * j2[1]; + + T d = j2[0] * j2[0] + j2[1] * j2[1]; + + vtkm::MatrixSetRow(jacobian, 0, vtkm::Vec(a, b)); + vtkm::MatrixSetRow(jacobian, 1, vtkm::Vec(b, d)); +} + +template +VTKM_EXEC_CONT void ComputeLeftCauchyGreenTensor(vtkm::Matrix& jacobian) +{ + vtkm::Vec j1 = vtkm::MatrixGetRow(jacobian, 0); + vtkm::Vec j2 = vtkm::MatrixGetRow(jacobian, 1); + vtkm::Vec j3 = vtkm::MatrixGetRow(jacobian, 2); + + // Left Cauchy Green Tensor is J*J^T + // j1[0] j1[1] j1[2] | j1[0] j2[0] j3[0] + // j2[0] j2[1] j2[2] | j1[1] j2[1] j3[1] + // j3[0] j3[1] j3[2] | j1[2] j2[2] j3[2] + + T a = j1[0] * j1[0] + j1[1] * j1[1] + j1[2] * j1[2]; + T b = j1[0] * j2[0] + j1[1] * j2[1] + j1[2] * j2[2]; + T c = j1[0] * j3[0] + j1[1] * j3[1] + j1[2] * j3[2]; + + T d = j2[0] * j2[0] + j2[1] * j2[1] + j2[2] * j2[2]; + T e = j2[0] * j3[0] + j2[1] * j3[1] + j2[2] * j3[2]; + + T f = j3[0] * j3[0] + j3[1] * j3[1] + j3[2] * j3[2]; + + vtkm::MatrixSetRow(jacobian, 0, vtkm::Vec(a, b, c)); + vtkm::MatrixSetRow(jacobian, 1, vtkm::Vec(b, d, e)); + vtkm::MatrixSetRow(jacobian, 2, vtkm::Vec(d, e, f)); +} + +template +VTKM_EXEC_CONT void ComputeRightCauchyGreenTensor(vtkm::Matrix& jacobian) +{ + vtkm::Vec j1 = vtkm::MatrixGetRow(jacobian, 0); + vtkm::Vec j2 = vtkm::MatrixGetRow(jacobian, 1); + + // Right Cauchy Green Tensor is J^T*J + // j1[0] j2[0] | j1[0] j1[1] + // j1[1] j2[1] | j2[0] j2[1] + + T a = j1[0] * j1[0] + j2[0] * j2[0]; + T b = j1[0] * j1[1] + j2[0] * j2[1]; + + T d = j1[1] * j1[1] + j2[1] * j2[1]; + + j1 = vtkm::Vec(a, b); + j2 = vtkm::Vec(b, d); +} + +template +VTKM_EXEC_CONT void ComputeRightCauchyGreenTensor(vtkm::Matrix& jacobian) +{ + vtkm::Vec j1 = vtkm::MatrixGetRow(jacobian, 0); + vtkm::Vec j2 = vtkm::MatrixGetRow(jacobian, 1); + vtkm::Vec j3 = vtkm::MatrixGetRow(jacobian, 2); + + // Right Cauchy Green Tensor is J^T*J + // j1[0] j2[0] j3[0] | j1[0] j1[1] j1[2] + // j1[1] j2[1] j3[1] | j2[0] j2[1] j2[2] + // j1[2] j2[2] j3[2] | j3[0] j3[1] j3[2] + + T a = j1[0] * j1[0] + j2[0] * j2[0] + j3[0] * j3[0]; + T b = j1[0] * j1[1] + j2[0] * j2[1] + j3[0] * j3[1]; + T c = j1[0] * j1[2] + j2[0] * j2[2] + j3[0] * j3[2]; + + T d = j1[1] * j1[1] + j2[1] * j2[1] + j3[1] * j3[1]; + T e = j1[1] * j1[2] + j2[1] * j2[2] + j3[1] * j3[2]; + + T f = j1[2] * j1[2] + j2[2] * j2[2] + j3[2] * j3[2]; + + j1 = vtkm::Vec(a, b, c); + j2 = vtkm::Vec(b, d, e); + j3 = vtkm::Vec(d, e, f); +} + +template +VTKM_EXEC_CONT void Jacobi(vtkm::Matrix tensor, vtkm::Vec& eigen) +{ + vtkm::Vec j1 = vtkm::MatrixGetRow(tensor, 0); + vtkm::Vec j2 = vtkm::MatrixGetRow(tensor, 1); + + // Assume a symetric matrix + // a b + // b c + + T a = j1[0]; + T b = j1[1]; + T c = j2[1]; + + T trace = (a + c) / 2.0f; + T det = a * c - b * b; + T sqrtr = vtkm::Sqrt(trace * trace - det); + + // Order the largest first to match VTK + eigen[0] = trace + sqrtr; + eigen[1] = trace - sqrtr; +} + +template +VTKM_EXEC_CONT inline void cswap(T& v1, T& v2) +{ + if (v2 < v1) + vtkm::Swap(v1, v2); +} + +template +VTKM_EXEC_CONT void Jacobi(vtkm::Matrix tensor, vtkm::Vec& eigen) +{ + vtkm::Vec j1 = vtkm::MatrixGetRow(tensor, 0); + vtkm::Vec j2 = vtkm::MatrixGetRow(tensor, 1); + vtkm::Vec j3 = vtkm::MatrixGetRow(tensor, 2); + + // Assume a symetric matrix + // a b c + // b d e + // c e f + T a = j1[0]; + T b = j1[1]; + T c = j1[2]; + T d = j2[1]; + T e = j2[2]; + T f = j3[2]; + + T x = (a + d + f) / 3.0f; // trace + + a -= x; + d -= x; + f -= x; + + // Det / 2; + T q = (a * d * f + b * e * c + c * b * e - c * d * c - e * e * a - f * b * b) / 2.0f; + T r = (a * a + b * b + c * c + b * b + d * d + e * e + c * c + e * e + f * f) / 6.0f; + + T D = (r * r * r - q * q); + T phi = 0.0f; + + if (D < vtkm::Epsilon()) + phi = 0.0f; + else + { + phi = vtkm::ATan(vtkm::Sqrt(D) / q) / 3.0f; + + if (phi < 0) + phi += static_cast(vtkm::Pi()); + } + + const T sqrt3 = vtkm::Sqrt(3.0f); + const T sqrtr = vtkm::Sqrt(r); + + T sinphi = 0.0f, cosphi = 0.0f; + sinphi = vtkm::Sin(phi); + cosphi = vtkm::Cos(phi); + + // Sorted in decreasing order. + T w0 = x + 2.0f * sqrtr * cosphi; + T w1 = x - sqrtr * (cosphi - sqrt3 * sinphi); + T w2 = x - sqrtr * (cosphi + sqrt3 * sinphi); + + cswap(w0, w1); + cswap(w0, w2); + cswap(w1, w2); + + eigen[0] = w0; + eigen[1] = w1; + eigen[2] = w2; +} + +} // namespace detail +} // namespace worklet +} // namespace vtkm +#endif //vtk_m_worklet_lcs_LagrangianStructureHelpers_h