diff --git a/.gitignore b/.gitignore new file mode 100644 index 000000000..e43b0f988 --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +.DS_Store diff --git a/docs/changelog/oscillator-source.md b/docs/changelog/oscillator-source.md new file mode 100644 index 000000000..e639e1601 --- /dev/null +++ b/docs/changelog/oscillator-source.md @@ -0,0 +1,11 @@ +# Time-varying "oscillator" source-filter and example + +The oscillator is a simple analytical source of time-varying data. +It provides a function value at each point of a uniform grid that +is computed as a sum of Gaussian kernels — each with a specified +position, amplitude, frequency, and phase. + +The example (in `examples/oscillator`) generates volumetric Cinema +datasets that can be viewed in a web browser with [ArcticViewer][]. + +[ArcticViewer]: https://kitware.github.io/arctic-viewer/ diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 4871f6e7d..f13e19d30 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -33,6 +33,7 @@ add_subdirectory(hello_world) add_subdirectory(histogram) add_subdirectory(isosurface) add_subdirectory(multi_backend) +add_subdirectory(oscillator) add_subdirectory(particle_advection) add_subdirectory(temporal_advection) add_subdirectory(redistribute_points) diff --git a/examples/oscillator/CMakeLists.txt b/examples/oscillator/CMakeLists.txt new file mode 100644 index 000000000..b3cb0d8d9 --- /dev/null +++ b/examples/oscillator/CMakeLists.txt @@ -0,0 +1,46 @@ +##============================================================================= +## +## Copyright (c) Kitware, Inc. +## All rights reserved. +## See LICENSE.txt for details. +## +## This software is distributed WITHOUT ANY WARRANTY; without even +## the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +## PURPOSE. See the above copyright notice for more information. +## +## Copyright 2015 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +## Copyright 2015 UT-Battelle, LLC. +## Copyright 2015 Los Alamos National Security. +## +## Under the terms of Contract DE-NA0003525 with NTESS, +## the U.S. Government retains certain rights in this software. +## Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National +## Laboratory (LANL), the U.S. Government retains certain rights in +## this software. +## +##============================================================================= +cmake_minimum_required(VERSION 3.3...3.12 FATAL_ERROR) +project(Oscillator CXX) + +#Find the VTK-m package +find_package(VTKm REQUIRED QUIET) + +add_executable(Oscillator_SERIAL Oscillator.cxx) +target_compile_definitions(Oscillator_SERIAL PRIVATE + "VTKM_DEVICE_ADAPTER=VTKM_DEVICE_ADAPTER_SERIAL") +target_link_libraries(Oscillator_SERIAL PRIVATE vtkm_cont) + +if (TARGET vtkm::tbb) + add_executable(Oscillator_TBB Oscillator.cxx) + target_compile_definitions(Oscillator_TBB PRIVATE + "VTKM_DEVICE_ADAPTER=VTKM_DEVICE_ADAPTER_TBB") + target_link_libraries(Oscillator_TBB PRIVATE vtkm_cont) +endif() + +if (TARGET vtkm::cuda) + vtkm_compile_as_cuda(oscillatorCudaSrc Oscillator.cxx) + add_executable(Oscillator_CUDA ${oscillatorCudaSrc}) + target_compile_definitions(Oscillator_CUDA PRIVATE + "VTKM_DEVICE_ADAPTER=VTKM_DEVICE_ADAPTER_CUDA") + target_link_libraries(Oscillator_CUDA PRIVATE vtkm_cont) +endif() diff --git a/examples/oscillator/Oscillator.cxx b/examples/oscillator/Oscillator.cxx new file mode 100644 index 000000000..412e96c58 --- /dev/null +++ b/examples/oscillator/Oscillator.cxx @@ -0,0 +1,357 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +// +// Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +// Copyright 2014 UT-Battelle, LLC. +// Copyright 2014 Los Alamos National Security. +// +// Under the terms of Contract DE-NA0003525 with NTESS, +// the U.S. Government retains certain rights in this software. +// +// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National +// Laboratory (LANL), the U.S. Government retains certain rights in +// this software. +//============================================================================ +#include +#include +#include +#include + +#include "vtkm/Math.h" +#include "vtkm/cont/ArrayHandle.h" +#include "vtkm/cont/DataSetBuilderUniform.h" + +#include "vtkm/filter/FilterDataSet.h" + +#include "vtkm/cont/TryExecute.h" +#include "vtkm/cont/cuda/DeviceAdapterCuda.h" +#include "vtkm/cont/serial/DeviceAdapterSerial.h" +#include "vtkm/cont/tbb/DeviceAdapterTBB.h" + +#include "vtkm/filter/OscillatorSource.h" + +#if !defined(_WIN32) || defined(__CYGWIN__) +#include /* unlink */ +#else +#include /* unlink */ +#endif + +//Suppress warnings about glut being deprecated on OSX +#if (defined(VTKM_GCC) || defined(VTKM_CLANG)) +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wdeprecated-declarations" +#endif + +//This is the list of devices to compile in support for. The order of the +//devices determines the runtime preference. +struct DevicesToTry : vtkm::ListTagBase +{ +}; + +// ---------------------------------------------------------------------------- + +struct OscillatorPolicy : public vtkm::filter::PolicyBase +{ + using DeviceAdapterList = DevicesToTry; +}; + +// trim() from http://stackoverflow.com/a/217605/44738 +static inline std::string& ltrim(std::string& s) +{ + s.erase(s.begin(), + std::find_if(s.begin(), s.end(), std::not1(std::ptr_fun(std::isspace)))); + return s; +} +static inline std::string& rtrim(std::string& s) +{ + s.erase( + std::find_if(s.rbegin(), s.rend(), std::not1(std::ptr_fun(std::isspace))).base(), + s.end()); + return s; +} +static inline std::string& trim(std::string& s) +{ + return ltrim(rtrim(s)); +} + +// ---------------------------------------------------------------------------- + +void read_oscillators(std::string filePath, vtkm::filter::OscillatorSource& filter) +{ + std::ifstream in(filePath); + if (!in) + throw std::runtime_error("Unable to open " + filePath); + std::string line; + while (std::getline(in, line)) + { + line = trim(line); + if (line.empty() || line[0] == '#') + continue; + std::istringstream iss(line); + + std::string stype; + iss >> stype; + + int x, y, z; + iss >> x >> y >> z; + + float r, omega0, zeta = 0.0f; + iss >> r >> omega0; + + if (stype == "damped") + { + iss >> zeta; + } + + if (stype == "damped") + { + filter.AddDamped(x, y, z, r, omega0, zeta); + } + else if (stype == "decaying") + { + filter.AddDecaying(x, y, z, r, omega0, zeta); + } + else if (stype == "periodic") + { + filter.AddPeriodic(x, y, z, r, omega0, zeta); + } + } +} + +// ---------------------------------------------------------------------------- +// ArcticViewer helper +// ---------------------------------------------------------------------------- + +void writeData(std::string& basePath, int timestep, int iSize, int jSize, int kSize, double* values) +{ + int size = iSize * jSize * kSize; + std::ostringstream timeValues; + for (int t = 0; t <= timestep; t++) + { + if (t > 0) + { + timeValues << ", "; + } + timeValues << "\"" << t << "\""; + } + + // Root metadata file + std::ostringstream metaFilePath; + metaFilePath << basePath.c_str() << "/index.json"; + std::ofstream metaFilePathPointer(metaFilePath.str().c_str(), std::ios::out); + if (metaFilePathPointer.fail()) + { + std::cout << "Unable to open file: " << metaFilePath.str().c_str() << std::endl; + } + else + { + metaFilePathPointer << "{\n" + << " \"metadata\": {\"backgroundColor\": \"#000000\"},\n" + << " \"type\": [\"tonic-query-data-model\", \"vtk-volume\"],\n" + << " \"arguments_order\": [\"time\"],\n" + << " \"arguments\": {\n" + << " \"time\": {\n" + << " \"loop\": \"modulo\",\n" + << " \"ui\": \"slider\",\n" + << " \"values\": [" << timeValues.str().c_str() << "]\n" + << " }\n" + << " },\n" + << " \"data\": [{\n" + << " \"pattern\": \"{time}.json\",\n" + << " \"rootFile\": true,\n" + << " \"name\": \"scene\",\n" + << " \"type\": \"json\"\n" + << " }]\n" + << "}\n"; + metaFilePathPointer.flush(); + metaFilePathPointer.close(); + } + + // Current data file + std::ostringstream dataFilePath; + dataFilePath << basePath.c_str() << "/" << timestep << ".data"; + std::ofstream dataFilePathPointer(dataFilePath.str().c_str(), std::ios::out | std::ios::binary); + if (dataFilePathPointer.fail()) + { + std::cout << "Unable to open file: " << dataFilePath.str().c_str() << std::endl; + } + else + { + int stackSize = size * 8; + dataFilePathPointer.write((char*)values, stackSize); + dataFilePathPointer.flush(); + dataFilePathPointer.close(); + } + + // Current dataset meta file + std::ostringstream dataMetaFilePath; + dataMetaFilePath << basePath.c_str() << "/" << timestep << ".json"; + std::ofstream dataMetaFilePathPointer(dataMetaFilePath.str().c_str(), std::ios::out); + if (dataMetaFilePathPointer.fail()) + { + std::cout << "Unable to open file: " << dataMetaFilePath.str().c_str() << std::endl; + } + else + { + dataMetaFilePathPointer << "{\n" + << " \"origin\": [0,0,0],\n" + << " \"spacing\": [1,1,1],\n" + << " \"extent\": [\n" + << " 0, " << (iSize - 1) << ",\n" + << " 0, " << (jSize - 1) << ",\n" + << " 0, " << (kSize - 1) << "],\n" + << " \"vtkClass\": \"vtkImageData\",\n" + << " \"pointData\": {\n" + << " \"vtkClass\": \"vtkDataSetAttributes\",\n" + << " \"arrays\": [{\n" + << " \"data\": {\n" + << " \"numberOfComponents\": 1,\n" + << " \"name\": \"oscillation\",\n" + << " \"vtkClass\": \"vtkDataArray\",\n" + << " \"dataType\": \"Float64Array\",\n" + << " \"ref\": {\n" + << " \"registration\": \"setScalars\",\n" + << " \"encode\": \"LittleEndian\",\n" + << " \"basepath\": \"\",\n" + << " \"id\": \"" << timestep << ".data\"\n" + << " },\n" + << " \"size\": " << size << std::endl + << " }\n" + << " }]\n" + << " }\n" + << "}\n"; + dataMetaFilePathPointer.flush(); + dataMetaFilePathPointer.close(); + } +} + +// ---------------------------------------------------------------------------- + +void printUsage() +{ + std::cout << "Usage: Oscillator [options]\n" + << std::endl + << "Options:\n" + << std::endl + << " -s, --shape POINT domain shape [default: 64x64x64]" << std::endl + << " -t, --dt FLOAT time step [default: 0.01]" << std::endl + << " -f, --config STRING oscillator file (required)" << std::endl + << " --t-end FLOAT end time [default: 10]" << std::endl + << " -o, --output STRING directory to output data\n" + << std::endl; +} + +// ---------------------------------------------------------------------------- + +int main(int argc, char** argv) +{ + std::string oscillatorConfigFile = ""; + std::string outputDirectory = ""; + int sizeX = 64; + int sizeY = 64; + int sizeZ = 64; + float startTime = 0.0f; + float endTime = 10.0f; + float deltaTime = 0.01f; + float currentTime = startTime; + bool generateOutput = false; + + // Process args + int nbOptions = argc - 1; + for (int i = 1; i < nbOptions; i += 2) + { + // Config (REQUIRED) + if (!strcmp(argv[i], "-f") || !strcmp(argv[i], "--config")) + { + oscillatorConfigFile = argv[i + 1]; + } + // Output + if (!strcmp(argv[i], "-o") || !strcmp(argv[i], "--output")) + { + generateOutput = true; + outputDirectory = argv[i + 1]; + } + + // Shape + if (!strcmp(argv[i], "-s") || !strcmp(argv[i], "--shape")) + { + std::string shape = argv[i + 1]; + size_t found = shape.find("x"); + sizeX = sizeY = sizeZ = std::stoi(shape); + if (found != std::string::npos) + { + shape = shape.substr(found + 1); + sizeY = sizeZ = std::stoi(shape); + found = shape.find("x"); + if (found != std::string::npos) + { + shape = shape.substr(found + 1); + sizeZ = std::stoi(shape); + found = shape.find("x"); + } + } + } + + // Time + if (!strcmp(argv[i], "-t") || !strcmp(argv[i], "--dt")) + { + deltaTime = float(std::atof(argv[i + 1])); + } + + // End-Time + if (!strcmp(argv[i], "--t-end")) + { + endTime = float(std::atof(argv[i + 1])); + } + } + + if (oscillatorConfigFile.size() < 2) + { + printUsage(); + return 0; + } + + std::cout << "\n=========== configuration ============" << std::endl; + std::cout << " oscillator config: " << oscillatorConfigFile.c_str() << std::endl; + std::cout << " mesh size: " << sizeX << "x" << sizeY << "x" << sizeZ << std::endl; + std::cout << " time handling:" << std::endl; + std::cout << " - dt: " << deltaTime << std::endl; + std::cout << " - end: " << endTime << std::endl; + std::cout << "=======================================\n" << std::endl; + + vtkm::cont::DataSetBuilderUniform builder; + vtkm::cont::DataSet dataset = builder.Create(vtkm::Id3(sizeX, sizeY, sizeZ)); + + vtkm::filter::OscillatorSource filter; + read_oscillators(oscillatorConfigFile, filter); + + std::cout << "=========== start computation ============" << std::endl; + int count = 0; + while (currentTime < endTime) + { + filter.SetTime(currentTime); + vtkm::cont::DataSet rdata = filter.Execute(dataset, OscillatorPolicy()); + if (generateOutput) + { + vtkm::cont::ArrayHandle tmp; + rdata.GetField("oscillation", vtkm::cont::Field::Association::POINTS).GetData().CopyTo(tmp); + double* values = tmp.GetStorage().GetArray(); + writeData(outputDirectory, count++, sizeX, sizeY, sizeZ, values); + } + + std::cout << "Compute time " << currentTime << std::endl; + currentTime += deltaTime; + } + std::cout << "=========== computation done ============" << std::endl; +} + +#if (defined(VTKM_GCC) || defined(VTKM_CLANG)) +#pragma GCC diagnostic pop +#endif diff --git a/examples/oscillator/inputs/complex-setup-256x256x128.osc b/examples/oscillator/inputs/complex-setup-256x256x128.osc new file mode 100644 index 000000000..142666f4f --- /dev/null +++ b/examples/oscillator/inputs/complex-setup-256x256x128.osc @@ -0,0 +1,14 @@ +# type center r omega0 zeta +periodic 44 44 64 35 0.785 +periodic 128 44 64 35 1.57 +periodic 212 44 64 35 2.355 +periodic 44 128 64 35 3.14 +periodic 128 128 64 35 3.925 +periodic 212 128 64 35 4.71 +periodic 44 212 64 35 5.495 +periodic 128 212 64 35 6.28 +periodic 212 212 64 35 7.065 +periodic 86 86 64 20 7.065 +periodic 86 170 64 20 5.495 +periodic 168 86 64 20 3.14 +periodic 170 170 64 20 1.57 diff --git a/examples/oscillator/inputs/periodic-64.osc b/examples/oscillator/inputs/periodic-64.osc new file mode 100644 index 000000000..ba8ce52e8 --- /dev/null +++ b/examples/oscillator/inputs/periodic-64.osc @@ -0,0 +1,9 @@ +# type center r omega0 zeta +periodic 16 16 16 8 3.14 +periodic 16 16 48 8 3.14 +periodic 16 48 16 8 3.14 +periodic 16 48 48 8 3.14 +periodic 48 16 16 8 3.14 +periodic 48 16 48 8 3.14 +periodic 48 48 16 8 3.14 +periodic 48 48 48 8 3.14 diff --git a/examples/oscillator/inputs/sample.osc b/examples/oscillator/inputs/sample.osc new file mode 100644 index 000000000..42650fc26 --- /dev/null +++ b/examples/oscillator/inputs/sample.osc @@ -0,0 +1,6 @@ +# type center r omega0 zeta +damped 32 32 32 10. 3.14 .3 +damped 16 32 16 10. 9.5 .1 +damped 48 32 48 5. 3.14 .1 +decaying 16 32 48 15 3.14 +periodic 48 32 16 15 3.14 diff --git a/vtkm/filter/CMakeLists.txt b/vtkm/filter/CMakeLists.txt index 69c5e38c8..c34691a24 100644 --- a/vtkm/filter/CMakeLists.txt +++ b/vtkm/filter/CMakeLists.txt @@ -49,6 +49,7 @@ set(headers MaskPoints.h NDEntropy.h NDHistogram.h + OscillatorSource.h PointAverage.h PointElevation.h PointTransform.h @@ -93,6 +94,7 @@ set(header_template_sources MaskPoints.hxx NDEntropy.hxx NDHistogram.hxx + OscillatorSource.hxx PointAverage.hxx PointElevation.hxx PointTransform.hxx diff --git a/vtkm/filter/OscillatorSource.h b/vtkm/filter/OscillatorSource.h new file mode 100644 index 000000000..fa671c2dd --- /dev/null +++ b/vtkm/filter/OscillatorSource.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. +// +// Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +// Copyright 2014 UT-Battelle, LLC. +// Copyright 2014 Los Alamos National Security. +// +// Under the terms of Contract DE-NA0003525 with NTESS, +// the U.S. Government retains certain rights in this software. +// +// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National +// Laboratory (LANL), the U.S. Government retains certain rights in +// this software. +//============================================================================ +#ifndef vtk_m_filter_OscillatorSource_h +#define vtk_m_filter_OscillatorSource_h + +#include "vtkm/filter/FilterField.h" +#include "vtkm/worklet/OscillatorSource.h" + +namespace vtkm +{ +namespace filter +{ + +/**\brief An analytical, time-varying array-source filter. + * + * This filter will create a new array (named "oscillation" by default) + * that evaluates to a sum of time-varying Gaussian exponentials + * specified in its configuration. + */ +class OscillatorSource : public vtkm::filter::FilterField +{ +public: + VTKM_CONT + OscillatorSource(); + + VTKM_CONT + void SetTime(vtkm::Float64 time); + + VTKM_CONT + void AddPeriodic(vtkm::Float64 x, + vtkm::Float64 y, + vtkm::Float64 z, + vtkm::Float64 radius, + vtkm::Float64 omega, + vtkm::Float64 zeta); + + VTKM_CONT + void AddDamped(vtkm::Float64 x, + vtkm::Float64 y, + vtkm::Float64 z, + vtkm::Float64 radius, + vtkm::Float64 omega, + vtkm::Float64 zeta); + + VTKM_CONT + void AddDecaying(vtkm::Float64 x, + vtkm::Float64 y, + vtkm::Float64 z, + vtkm::Float64 radius, + vtkm::Float64 omega, + vtkm::Float64 zeta); + + template + VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input, + const vtkm::cont::ArrayHandle& field, + const vtkm::filter::FieldMetadata& fieldMeta, + const vtkm::filter::PolicyBase& policy, + const DeviceAdapter& tag); + +private: + vtkm::worklet::OscillatorSource Worklet; +}; + +template <> +class FilterTraits +{ +public: + //Point Oscillator can only convert Float and Double Vec3 arrays + using InputFieldTypeList = vtkm::TypeListTagFieldVec3; +}; +} +} + +#include "vtkm/filter/OscillatorSource.hxx" + +#endif // vtk_m_filter_OscillatorSource_h diff --git a/vtkm/filter/OscillatorSource.hxx b/vtkm/filter/OscillatorSource.hxx new file mode 100644 index 000000000..4c035fdf9 --- /dev/null +++ b/vtkm/filter/OscillatorSource.hxx @@ -0,0 +1,101 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +// +// Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +// Copyright 2014 UT-Battelle, LLC. +// Copyright 2014 Los Alamos National Security. +// +// Under the terms of Contract DE-NA0003525 with NTESS, +// the U.S. Government retains certain rights in this software. +// +// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National +// Laboratory (LANL), the U.S. Government retains certain rights in +// this software. +//============================================================================ +#include + +#include + +namespace vtkm +{ +namespace filter +{ + +//----------------------------------------------------------------------------- +inline VTKM_CONT OscillatorSource::OscillatorSource() + : Worklet() +{ + this->SetUseCoordinateSystemAsField(true); + this->SetOutputFieldName("oscillation"); +} + +//----------------------------------------------------------------------------- +inline VTKM_CONT void OscillatorSource::SetTime(vtkm::Float64 time) +{ + this->Worklet.SetTime(time); +} + +//----------------------------------------------------------------------------- +inline VTKM_CONT void OscillatorSource::AddPeriodic(vtkm::Float64 x, + vtkm::Float64 y, + vtkm::Float64 z, + vtkm::Float64 radius, + vtkm::Float64 omega, + vtkm::Float64 zeta) +{ + this->Worklet.AddPeriodic(x, y, z, radius, omega, zeta); +} + +//----------------------------------------------------------------------------- +inline VTKM_CONT void OscillatorSource::AddDamped(vtkm::Float64 x, + vtkm::Float64 y, + vtkm::Float64 z, + vtkm::Float64 radius, + vtkm::Float64 omega, + vtkm::Float64 zeta) +{ + this->Worklet.AddDamped(x, y, z, radius, omega, zeta); +} + +//----------------------------------------------------------------------------- +inline VTKM_CONT void OscillatorSource::AddDecaying(vtkm::Float64 x, + vtkm::Float64 y, + vtkm::Float64 z, + vtkm::Float64 radius, + vtkm::Float64 omega, + vtkm::Float64 zeta) +{ + this->Worklet.AddDecaying(x, y, z, radius, omega, zeta); +} + + +//----------------------------------------------------------------------------- +template +inline VTKM_CONT vtkm::cont::DataSet OscillatorSource::DoExecute( + const vtkm::cont::DataSet& inDataSet, + const vtkm::cont::ArrayHandle& field, + const vtkm::filter::FieldMetadata& fieldMetadata, + const vtkm::filter::PolicyBase&, + const DeviceAdapter&) +{ + vtkm::cont::ArrayHandle outArray; + vtkm::worklet::DispatcherMapField dispatcher( + this->Worklet); + + //todo, we need to use the policy to determine the valid conversions + //that the dispatcher should do + dispatcher.Invoke(field, outArray); + + return internal::CreateResult(inDataSet, + outArray, + this->GetOutputFieldName(), + fieldMetadata.GetAssociation(), + fieldMetadata.GetCellSetName()); +} +} +} // namespace vtkm::filter diff --git a/vtkm/worklet/CMakeLists.txt b/vtkm/worklet/CMakeLists.txt index ac3120821..07caab626 100644 --- a/vtkm/worklet/CMakeLists.txt +++ b/vtkm/worklet/CMakeLists.txt @@ -53,6 +53,7 @@ set(headers NDimsHistogram.h NDimsHistMarginalization.h Normalize.h + OscillatorSource.h ParticleAdvection.h PointAverage.h PointElevation.h diff --git a/vtkm/worklet/OscillatorSource.h b/vtkm/worklet/OscillatorSource.h new file mode 100644 index 000000000..6d3d6b5d8 --- /dev/null +++ b/vtkm/worklet/OscillatorSource.h @@ -0,0 +1,197 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +// +// Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +// Copyright 2014 UT-Battelle, LLC. +// Copyright 2014 Los Alamos National Security. +// +// Under the terms of Contract DE-NA0003525 with NTESS, +// the U.S. Government retains certain rights in this software. +// +// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National +// Laboratory (LANL), the U.S. Government retains certain rights in +// this software. +//============================================================================ +#ifndef vtk_m_worklet_OscillatorSource_h +#define vtk_m_worklet_OscillatorSource_h + +#include +#include + +#define MAX_OSCILLATORS 10 + +namespace vtkm +{ +namespace worklet +{ +namespace internal +{ + +struct Oscillator +{ + void Set(vtkm::Float64 x, + vtkm::Float64 y, + vtkm::Float64 z, + vtkm::Float64 radius, + vtkm::Float64 omega, + vtkm::Float64 zeta) + { + this->Center[0] = x; + this->Center[1] = y; + this->Center[2] = z; + this->Radius = radius; + this->Omega = omega; + this->Zeta = zeta; + } + + vtkm::Vec Center; + vtkm::Float64 Radius; + vtkm::Float64 Omega; + vtkm::Float64 Zeta; +}; +} + +class OscillatorSource : public vtkm::worklet::WorkletMapField +{ +public: + typedef void ControlSignature(FieldIn, FieldOut); + typedef _2 ExecutionSignature(_1); + + VTKM_CONT + OscillatorSource() + : NumberOfPeriodics(0) + , NumberOfDamped(0) + , NumberOfDecaying(0) + { + } + + VTKM_CONT + void AddPeriodic(vtkm::Float64 x, + vtkm::Float64 y, + vtkm::Float64 z, + vtkm::Float64 radius, + vtkm::Float64 omega, + vtkm::Float64 zeta) + { + if (this->NumberOfPeriodics < MAX_OSCILLATORS) + { + this->PeriodicOscillators[this->NumberOfPeriodics].Set(x, y, z, radius, omega, zeta); + this->NumberOfPeriodics++; + } + } + + VTKM_CONT + void AddDamped(vtkm::Float64 x, + vtkm::Float64 y, + vtkm::Float64 z, + vtkm::Float64 radius, + vtkm::Float64 omega, + vtkm::Float64 zeta) + { + if (this->NumberOfDamped < MAX_OSCILLATORS) + { + this->DampedOscillators[this->NumberOfDamped * 6].Set(x, y, z, radius, omega, zeta); + this->NumberOfDamped++; + } + } + + VTKM_CONT + void AddDecaying(vtkm::Float64 x, + vtkm::Float64 y, + vtkm::Float64 z, + vtkm::Float64 radius, + vtkm::Float64 omega, + vtkm::Float64 zeta) + { + if (this->NumberOfDecaying < MAX_OSCILLATORS) + { + this->DecayingOscillators[this->NumberOfDecaying * 6].Set(x, y, z, radius, omega, zeta); + this->NumberOfDecaying++; + } + } + + VTKM_CONT + void SetTime(vtkm::Float64 time) { this->Time = time; } + + VTKM_EXEC + vtkm::Float64 operator()(const vtkm::Vec& vec) const + { + vtkm::UInt8 oIdx; + vtkm::Float64 t0, t, result = 0; + const internal::Oscillator* oscillator; + + t0 = 0.0; + t = this->Time * 2 * 3.14159265358979323846; + + // Compute damped + for (oIdx = 0; oIdx < this->NumberOfDamped; oIdx++) + { + oscillator = &this->DampedOscillators[oIdx]; + + vtkm::Vec delta = oscillator->Center - vec; + vtkm::Float64 dist2 = dot(delta, delta); + vtkm::Float64 dist_damp = vtkm::Exp(-dist2 / (2 * oscillator->Radius * oscillator->Radius)); + vtkm::Float64 phi = vtkm::ACos(oscillator->Zeta); + vtkm::Float64 val = 1. - + vtkm::Exp(-oscillator->Zeta * oscillator->Omega * t0) * + (vtkm::Sin(vtkm::Sqrt(1 - oscillator->Zeta * oscillator->Zeta) * oscillator->Omega * t + + phi) / + vtkm::Sin(phi)); + result += val * dist_damp; + } + + // Compute decaying + for (oIdx = 0; oIdx < this->NumberOfDecaying; oIdx++) + { + oscillator = &this->DecayingOscillators[oIdx]; + t = t0 + 1 / oscillator->Omega; + vtkm::Vec delta = oscillator->Center - vec; + vtkm::Float64 dist2 = dot(delta, delta); + vtkm::Float64 dist_damp = vtkm::Exp(-dist2 / (2 * oscillator->Radius * oscillator->Radius)); + vtkm::Float64 val = vtkm::Sin(t / oscillator->Omega) / (oscillator->Omega * t); + result += val * dist_damp; + } + + // Compute periodic + for (oIdx = 0; oIdx < this->NumberOfPeriodics; oIdx++) + { + oscillator = &this->PeriodicOscillators[oIdx]; + t = t0 + 1 / oscillator->Omega; + vtkm::Vec delta = oscillator->Center - vec; + vtkm::Float64 dist2 = dot(delta, delta); + vtkm::Float64 dist_damp = vtkm::Exp(-dist2 / (2 * oscillator->Radius * oscillator->Radius)); + vtkm::Float64 val = vtkm::Sin(t / oscillator->Omega); + result += val * dist_damp; + } + + // We are done... + return result; + } + + template + VTKM_EXEC vtkm::Float64 operator()(const vtkm::Vec& vec) const + { + return (*this)(vtkm::make_Vec(static_cast(vec[0]), + static_cast(vec[1]), + static_cast(vec[2]))); + } + +private: + vtkm::Vec PeriodicOscillators; + vtkm::Vec DampedOscillators; + vtkm::Vec DecayingOscillators; + vtkm::UInt8 NumberOfPeriodics; + vtkm::UInt8 NumberOfDamped; + vtkm::UInt8 NumberOfDecaying; + vtkm::Float64 Time; +}; +} + +} // namespace vtkm + +#endif // vtk_m_worklet_PointElevation_h