Merge topic 'lagrangian-coherent-structures'
ac8d8c95b Separating GridMetaData class 4aa51fcab Adding custom dataset name to LCS filter output ec8a29e07 Resolving nvcc implicit conversion warning 469e60430 Resolving nvcc compiler warnings d5ef47040 LCS fixes da3983d21 Changes from Ken's review 6b8aa3f24 Fixing missing header in install issue 9634924e1 Fixing issues with LCS Filter files ... Acked-by: Kitware Robot <kwrobot@kitware.com> Merge-request: !1764
This commit is contained in:
Abhishek Yenpure
Kitware Robot
commit
46d4e50622
13
docs/changelog/lagrangian-coherent-structures.md
Normal file
13
docs/changelog/lagrangian-coherent-structures.md
Normal file
@ -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.
|
||||
25
examples/lagrangian_structures/CMakeLists.txt
Normal file
25
examples/lagrangian_structures/CMakeLists.txt
Normal file
@ -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()
|
||||
58
examples/lagrangian_structures/LagrangianStructures.cxx
Normal file
58
examples/lagrangian_structures/LagrangianStructures.cxx
Normal file
@ -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 <cmath>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
#include <vtkm/cont/ArrayCopy.h>
|
||||
#include <vtkm/cont/DataSet.h>
|
||||
#include <vtkm/cont/DataSetBuilderUniform.h>
|
||||
#include <vtkm/cont/DataSetFieldAdd.h>
|
||||
#include <vtkm/cont/Initialize.h>
|
||||
|
||||
#include <vtkm/io/reader/VTKDataSetReader.h>
|
||||
#include <vtkm/io/writer/VTKDataSetWriter.h>
|
||||
|
||||
#include <vtkm/filter/LagrangianStructures.h>
|
||||
|
||||
int main(int argc, char** argv)
|
||||
{
|
||||
vtkm::cont::Initialize(argc, argv);
|
||||
|
||||
if (argc < 3)
|
||||
{
|
||||
std::cout << "Usage : flte <input dataset> <vector field name>" << 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;
|
||||
}
|
||||
@ -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
|
||||
|
||||
92
vtkm/filter/LagrangianStructures.h
Normal file
92
vtkm/filter/LagrangianStructures.h
Normal file
@ -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 <vtkm/filter/FilterDataSetWithField.h>
|
||||
#include <vtkm/worklet/ParticleAdvection.h>
|
||||
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
|
||||
#include <vtkm/worklet/particleadvection/Integrators.h>
|
||||
|
||||
namespace vtkm
|
||||
{
|
||||
namespace filter
|
||||
{
|
||||
|
||||
class LagrangianStructures : public vtkm::filter::FilterDataSetWithField<LagrangianStructures>
|
||||
{
|
||||
public:
|
||||
using SupportedTypes = vtkm::TypeListTagFieldVec3;
|
||||
|
||||
using Scalar = vtkm::worklet::particleadvection::ScalarType;
|
||||
using Vector = vtkm::Vec<Scalar, 3>;
|
||||
|
||||
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<Vector>& flowMap)
|
||||
{
|
||||
this->FlowMapOutput = flowMap;
|
||||
}
|
||||
inline vtkm::cont::ArrayHandle<Vector> GetFlowMapOutput() { return this->FlowMapOutput; }
|
||||
|
||||
template <typename T, typename StorageType, typename DerivedPolicy>
|
||||
VTKM_CONT vtkm::cont::DataSet DoExecute(
|
||||
const vtkm::cont::DataSet& input,
|
||||
const vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, StorageType>& field,
|
||||
const vtkm::filter::FieldMetadata& fieldMeta,
|
||||
const vtkm::filter::PolicyBase<DerivedPolicy>& policy);
|
||||
|
||||
//Map a new field onto the resulting dataset after running the filter
|
||||
//this call is only valid
|
||||
template <typename T, typename StorageType, typename DerivedPolicy>
|
||||
VTKM_CONT bool DoMapField(vtkm::cont::DataSet& result,
|
||||
const vtkm::cont::ArrayHandle<T, StorageType>& input,
|
||||
const vtkm::filter::FieldMetadata& fieldMeta,
|
||||
vtkm::filter::PolicyBase<DerivedPolicy> policy);
|
||||
|
||||
private:
|
||||
Scalar StepSize;
|
||||
vtkm::Id NumberOfSteps;
|
||||
Scalar AdvectionTime;
|
||||
bool UseAuxiliaryGrid = false;
|
||||
vtkm::Id3 AuxiliaryDims;
|
||||
bool UseFlowMapOutput = false;
|
||||
std::string OutputFieldName;
|
||||
vtkm::cont::ArrayHandle<Vector> FlowMapOutput;
|
||||
};
|
||||
|
||||
} // namespace filter
|
||||
} // namespace vtkm
|
||||
|
||||
#include <vtkm/filter/LagrangianStructures.hxx>
|
||||
|
||||
#endif // vtk_m_filter_LagrangianStructures_h
|
||||
146
vtkm/filter/LagrangianStructures.hxx
Normal file
146
vtkm/filter/LagrangianStructures.hxx
Normal file
@ -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 <vtkm/cont/ArrayCopy.h>
|
||||
#include <vtkm/cont/ArrayHandleIndex.h>
|
||||
#include <vtkm/cont/ErrorFilterExecution.h>
|
||||
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
|
||||
#include <vtkm/worklet/particleadvection/Integrators.h>
|
||||
#include <vtkm/worklet/particleadvection/Particles.h>
|
||||
|
||||
#include <vtkm/worklet/LagrangianStructures.h>
|
||||
|
||||
namespace vtkm
|
||||
{
|
||||
namespace filter
|
||||
{
|
||||
|
||||
//-----------------------------------------------------------------------------
|
||||
inline VTKM_CONT LagrangianStructures::LagrangianStructures()
|
||||
: vtkm::filter::FilterDataSetWithField<LagrangianStructures>()
|
||||
{
|
||||
OutputFieldName = std::string("FTLE");
|
||||
}
|
||||
|
||||
//-----------------------------------------------------------------------------
|
||||
template <typename T, typename StorageType, typename DerivedPolicy>
|
||||
inline VTKM_CONT vtkm::cont::DataSet LagrangianStructures::DoExecute(
|
||||
const vtkm::cont::DataSet& input,
|
||||
const vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, StorageType>& field,
|
||||
const vtkm::filter::FieldMetadata& fieldMeta,
|
||||
const vtkm::filter::PolicyBase<DerivedPolicy>&)
|
||||
{
|
||||
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<vtkm::Vec<T, 3>, StorageType>;
|
||||
using GridEvaluator = vtkm::worklet::particleadvection::GridEvaluator<FieldHandle>;
|
||||
using Integrator = vtkm::worklet::particleadvection::RK4Integrator<GridEvaluator>;
|
||||
|
||||
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<Scalar>(inputBounds.X.Min),
|
||||
static_cast<Scalar>(inputBounds.Y.Min),
|
||||
static_cast<Scalar>(inputBounds.Z.Min));
|
||||
Vector spacing;
|
||||
spacing[0] =
|
||||
static_cast<Scalar>(inputBounds.X.Length()) / static_cast<Scalar>(lcsGridDims[0] - 1);
|
||||
spacing[1] =
|
||||
static_cast<Scalar>(inputBounds.Y.Length()) / static_cast<Scalar>(lcsGridDims[1] - 1);
|
||||
spacing[2] =
|
||||
static_cast<Scalar>(inputBounds.Z.Length()) / static_cast<Scalar>(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<Structured2DType>() || cellset.IsType<Structured3DType>()))
|
||||
throw vtkm::cont::ErrorFilterExecution(
|
||||
"Provided data is not structured, provide parameters for an auxiliary grid.");
|
||||
lcsInput = input;
|
||||
}
|
||||
vtkm::cont::ArrayHandle<Vector> 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<Vector> advectionPoints;
|
||||
vtkm::cont::ArrayCopy(lcsInputPoints, advectionPoints);
|
||||
advectionResult = particles.Run(integrator, advectionPoints, numberOfSteps);
|
||||
lcsOutputPoints = advectionResult.positions;
|
||||
}
|
||||
// FTLE output field
|
||||
vtkm::cont::ArrayHandle<Scalar> outputField;
|
||||
vtkm::FloatDefault advectionTime = this->GetAdvectionTime();
|
||||
|
||||
vtkm::cont::DynamicCellSet lcsCellSet = lcsInput.GetCellSet();
|
||||
if (lcsCellSet.IsType<Structured2DType>())
|
||||
{
|
||||
using AnalysisType = vtkm::worklet::LagrangianStructures<2>;
|
||||
AnalysisType ftleCalculator(advectionTime, lcsCellSet);
|
||||
vtkm::worklet::DispatcherMapField<AnalysisType> dispatcher(ftleCalculator);
|
||||
dispatcher.Invoke(lcsInputPoints, lcsOutputPoints, outputField);
|
||||
}
|
||||
else if (lcsCellSet.IsType<Structured3DType>())
|
||||
{
|
||||
using AnalysisType = vtkm::worklet::LagrangianStructures<3>;
|
||||
AnalysisType ftleCalculator(advectionTime, lcsCellSet);
|
||||
vtkm::worklet::DispatcherMapField<AnalysisType> 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 <typename T, typename StorageType, typename DerivedPolicy>
|
||||
inline VTKM_CONT bool LagrangianStructures::DoMapField(
|
||||
vtkm::cont::DataSet&,
|
||||
const vtkm::cont::ArrayHandle<T, StorageType>&,
|
||||
const vtkm::filter::FieldMetadata&,
|
||||
vtkm::filter::PolicyBase<DerivedPolicy>)
|
||||
{
|
||||
return false;
|
||||
}
|
||||
}
|
||||
} // namespace vtkm::filter
|
||||
@ -34,6 +34,7 @@ set(unit_tests
|
||||
UnitTestHistogramFilter.cxx
|
||||
UnitTestImageConnectivityFilter.cxx
|
||||
UnitTestLagrangianFilter.cxx
|
||||
UnitTestLagrangianStructuresFilter.cxx
|
||||
UnitTestMarchingCubesFilter.cxx
|
||||
UnitTestMaskFilter.cxx
|
||||
UnitTestMaskPointsFilter.cxx
|
||||
|
||||
228
vtkm/filter/testing/UnitTestLagrangianStructuresFilter.cxx
Normal file
228
vtkm/filter/testing/UnitTestLagrangianStructuresFilter.cxx
Normal file
@ -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 <iostream>
|
||||
#include <vtkm/cont/DataSetBuilderUniform.h>
|
||||
#include <vtkm/cont/DataSetFieldAdd.h>
|
||||
#include <vtkm/cont/testing/Testing.h>
|
||||
|
||||
#include <vtkm/filter/LagrangianStructures.h>
|
||||
|
||||
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<Scalar, 3>;
|
||||
|
||||
void ValidateLCSFilterResult(const vtkm::cont::ArrayHandle<Scalar>& vtkmOutput,
|
||||
const vtkm::cont::ArrayHandle<Scalar>& visitOutput,
|
||||
const vtkm::cont::ArrayHandle<Scalar>& 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<Scalar>()),
|
||||
"Calculated o/p is not as expected");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void TestLagrangianStructures()
|
||||
{
|
||||
using Scalar = vtkm::FloatDefault;
|
||||
using Field = vtkm::Vec<Scalar, 3>;
|
||||
|
||||
/*Construct dataset and vector field*/
|
||||
vtkm::cont::DataSet inputData;
|
||||
vtkm::Id2 dims(20, 10);
|
||||
vtkm::Vec<Scalar, 2> origin(0.0f, 0.0f);
|
||||
vtkm::Vec<Scalar, 2> spacing;
|
||||
spacing[0] = 2.0f / static_cast<Scalar>(dims[0] - 1);
|
||||
spacing[1] = 1.0f / static_cast<Scalar>(dims[1] - 1);
|
||||
vtkm::cont::DataSetBuilderUniform dataBuilder;
|
||||
inputData = dataBuilder.Create(dims, origin, spacing);
|
||||
|
||||
std::vector<Scalar> diffVec;
|
||||
std::vector<Scalar> visitVec;
|
||||
std::vector<Field> 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<Scalar> diffHandle = vtkm::cont::make_ArrayHandle(diffVec);
|
||||
vtkm::cont::ArrayHandle<Scalar> visitHandle = vtkm::cont::make_ArrayHandle(visitVec);
|
||||
vtkm::cont::ArrayHandle<Field> 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<Scalar> 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);
|
||||
}
|
||||
@ -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)
|
||||
|
||||
196
vtkm/worklet/LagrangianStructures.h
Normal file
196
vtkm/worklet/LagrangianStructures.h
Normal file
@ -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 <vtkm/Matrix.h>
|
||||
#include <vtkm/Types.h>
|
||||
#include <vtkm/worklet/DispatcherMapField.h>
|
||||
#include <vtkm/worklet/WorkletMapField.h>
|
||||
|
||||
#include <vtkm/worklet/lcs/GridMetaData.h>
|
||||
#include <vtkm/worklet/lcs/LagrangianStructureHelpers.h>
|
||||
|
||||
namespace vtkm
|
||||
{
|
||||
namespace worklet
|
||||
{
|
||||
|
||||
template <vtkm::IdComponent dimensions>
|
||||
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 <typename PointArray>
|
||||
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<vtkm::Id, 6> 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<Scalar, 2, 2> jacobian;
|
||||
vtkm::MatrixSetRow(jacobian, 0, vtkm::Vec<Scalar, 2>(f1x, f1y));
|
||||
vtkm::MatrixSetRow(jacobian, 1, vtkm::Vec<Scalar, 2>(f2x, f2y));
|
||||
|
||||
detail::ComputeLeftCauchyGreenTensor(jacobian);
|
||||
|
||||
vtkm::Vec<Scalar, 2> 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<Scalar>(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 <typename PointArray>
|
||||
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<vtkm::Id, 6> 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<Scalar, 3, 3> jacobian;
|
||||
vtkm::MatrixSetRow(jacobian, 0, vtkm::Vec<Scalar, 3>(f1x, f1y, f1z));
|
||||
vtkm::MatrixSetRow(jacobian, 1, vtkm::Vec<Scalar, 3>(f2x, f2y, f2z));
|
||||
vtkm::MatrixSetRow(jacobian, 2, vtkm::Vec<Scalar, 3>(f3x, f3y, f3z));
|
||||
|
||||
detail::ComputeLeftCauchyGreenTensor(jacobian);
|
||||
|
||||
vtkm::Vec<Scalar, 3> 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<Scalar>(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
|
||||
16
vtkm/worklet/lcs/CMakeLists.txt
Normal file
16
vtkm/worklet/lcs/CMakeLists.txt
Normal file
@ -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})
|
||||
93
vtkm/worklet/lcs/GridMetaData.h
Normal file
93
vtkm/worklet/lcs/GridMetaData.h
Normal file
@ -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<Structured2DType>())
|
||||
{
|
||||
this->cellSet2D = true;
|
||||
vtkm::Id2 dims =
|
||||
cellSet.Cast<Structured2DType>().GetSchedulingRange(vtkm::TopologyElementTagPoint());
|
||||
this->Dims = vtkm::Id3(dims[0], dims[1], 1);
|
||||
}
|
||||
else
|
||||
{
|
||||
this->cellSet2D = false;
|
||||
this->Dims =
|
||||
cellSet.Cast<Structured3DType>().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<vtkm::Id, 6> GetNeighborIndices(const vtkm::Id index) const
|
||||
{
|
||||
vtkm::Vec<vtkm::Id, 6> 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
|
||||
209
vtkm/worklet/lcs/LagrangianStructureHelpers.h
Normal file
209
vtkm/worklet/lcs/LagrangianStructureHelpers.h
Normal file
@ -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 <vtkm/Matrix.h>
|
||||
#include <vtkm/Swap.h>
|
||||
#include <vtkm/Types.h>
|
||||
|
||||
namespace vtkm
|
||||
{
|
||||
namespace worklet
|
||||
{
|
||||
namespace detail
|
||||
{
|
||||
|
||||
template <typename T>
|
||||
VTKM_EXEC_CONT void ComputeLeftCauchyGreenTensor(vtkm::Matrix<T, 2, 2>& jacobian)
|
||||
{
|
||||
vtkm::Vec<T, 2> j1 = vtkm::MatrixGetRow(jacobian, 0);
|
||||
vtkm::Vec<T, 2> 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<T, 2>(a, b));
|
||||
vtkm::MatrixSetRow(jacobian, 1, vtkm::Vec<T, 2>(b, d));
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
VTKM_EXEC_CONT void ComputeLeftCauchyGreenTensor(vtkm::Matrix<T, 3, 3>& jacobian)
|
||||
{
|
||||
vtkm::Vec<T, 3> j1 = vtkm::MatrixGetRow(jacobian, 0);
|
||||
vtkm::Vec<T, 3> j2 = vtkm::MatrixGetRow(jacobian, 1);
|
||||
vtkm::Vec<T, 3> 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<T, 3>(a, b, c));
|
||||
vtkm::MatrixSetRow(jacobian, 1, vtkm::Vec<T, 3>(b, d, e));
|
||||
vtkm::MatrixSetRow(jacobian, 2, vtkm::Vec<T, 3>(d, e, f));
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
VTKM_EXEC_CONT void ComputeRightCauchyGreenTensor(vtkm::Matrix<T, 2, 2>& jacobian)
|
||||
{
|
||||
vtkm::Vec<T, 2> j1 = vtkm::MatrixGetRow(jacobian, 0);
|
||||
vtkm::Vec<T, 2> 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<T, 2>(a, b);
|
||||
j2 = vtkm::Vec<T, 2>(b, d);
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
VTKM_EXEC_CONT void ComputeRightCauchyGreenTensor(vtkm::Matrix<T, 3, 3>& jacobian)
|
||||
{
|
||||
vtkm::Vec<T, 3> j1 = vtkm::MatrixGetRow(jacobian, 0);
|
||||
vtkm::Vec<T, 3> j2 = vtkm::MatrixGetRow(jacobian, 1);
|
||||
vtkm::Vec<T, 3> 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<T, 3>(a, b, c);
|
||||
j2 = vtkm::Vec<T, 3>(b, d, e);
|
||||
j3 = vtkm::Vec<T, 3>(d, e, f);
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
VTKM_EXEC_CONT void Jacobi(vtkm::Matrix<T, 2, 2> tensor, vtkm::Vec<T, 2>& eigen)
|
||||
{
|
||||
vtkm::Vec<T, 2> j1 = vtkm::MatrixGetRow(tensor, 0);
|
||||
vtkm::Vec<T, 2> 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 <typename T>
|
||||
VTKM_EXEC_CONT inline void cswap(T& v1, T& v2)
|
||||
{
|
||||
if (v2 < v1)
|
||||
vtkm::Swap(v1, v2);
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
VTKM_EXEC_CONT void Jacobi(vtkm::Matrix<T, 3, 3> tensor, vtkm::Vec<T, 3>& eigen)
|
||||
{
|
||||
vtkm::Vec<T, 3> j1 = vtkm::MatrixGetRow(tensor, 0);
|
||||
vtkm::Vec<T, 3> j2 = vtkm::MatrixGetRow(tensor, 1);
|
||||
vtkm::Vec<T, 3> 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<T>())
|
||||
phi = 0.0f;
|
||||
else
|
||||
{
|
||||
phi = vtkm::ATan(vtkm::Sqrt(D) / q) / 3.0f;
|
||||
|
||||
if (phi < 0)
|
||||
phi += static_cast<T>(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
|
||||
Loading…
Reference in New Issue
Block a user