Shuenhoy/vtk-m2
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add a point-oscillator filter + example

Browse Source 
The oscillator is a simple analytical source of time-varying data.
It provides a function value at each point that is computed as a
sum of Gaussian kernels -- each with a specified position, amplitude,
frequency, and phase.
This commit is contained in:
Sebastien Jourdain 2017-08-30 14:27:04 -06:00 committed by David Thompson
parent 776d085464
commit 4192b9a1d8
13 changed files with 839 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitignore vendored Normal file
View File

@ -0,0 +1 @@
.DS_Store

11
docs/changelog/oscillator-source.md Normal file
View File

@ -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/

1
examples/CMakeLists.txt
View File

@ -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)

46
examples/oscillator/CMakeLists.txt Normal file
View File

@ -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()

357
examples/oscillator/Oscillator.cxx Normal file
View File

@ -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 <algorithm>
#include <fstream>
#include <iostream>
#include <sstream>
#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 <unistd.h> /* unlink */
#else
#include <io.h> /* 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<vtkm::cont::DeviceAdapterTagCuda,
vtkm::cont::DeviceAdapterTagTBB,
vtkm::cont::DeviceAdapterTagSerial>
{
};
// ----------------------------------------------------------------------------
struct OscillatorPolicy : public vtkm::filter::PolicyBase<OscillatorPolicy>
{
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<int, int>(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<int, int>(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<vtkm::Float64> 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

14
examples/oscillator/inputs/complex-setup-256x256x128.osc Normal file
View File

@ -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

9
examples/oscillator/inputs/periodic-64.osc Normal file
View File

@ -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

6
examples/oscillator/inputs/sample.osc Normal file
View File

@ -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

2
vtkm/filter/CMakeLists.txt
View File

@ -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

93
vtkm/filter/OscillatorSource.h Normal file
View 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.
//
// 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<OscillatorSource>
{
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 <typename T, typename StorageType, typename DerivedPolicy, typename DeviceAdapter>
VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input,
const vtkm::cont::ArrayHandle<T, StorageType>& field,
const vtkm::filter::FieldMetadata& fieldMeta,
const vtkm::filter::PolicyBase<DerivedPolicy>& policy,
const DeviceAdapter& tag);
private:
vtkm::worklet::OscillatorSource Worklet;
};
template <>
class FilterTraits<OscillatorSource>
{
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

101
vtkm/filter/OscillatorSource.hxx Normal file
View File

@ -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 <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/filter/internal/CreateResult.h>
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 <typename T, typename StorageType, typename DerivedPolicy, typename DeviceAdapter>
inline VTKM_CONT vtkm::cont::DataSet OscillatorSource::DoExecute(
const vtkm::cont::DataSet& inDataSet,
const vtkm::cont::ArrayHandle<T, StorageType>& field,
const vtkm::filter::FieldMetadata& fieldMetadata,
const vtkm::filter::PolicyBase<DerivedPolicy>&,
const DeviceAdapter&)
{
vtkm::cont::ArrayHandle<vtkm::Float64> outArray;
vtkm::worklet::DispatcherMapField<vtkm::worklet::OscillatorSource, DeviceAdapter> 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

1
vtkm/worklet/CMakeLists.txt
View File

@ -53,6 +53,7 @@ set(headers
NDimsHistogram.h
NDimsHistMarginalization.h
Normalize.h
OscillatorSource.h
ParticleAdvection.h
PointAverage.h
PointElevation.h

197
vtkm/worklet/OscillatorSource.h Normal file
View File

@ -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 <vtkm/Math.h>
#include <vtkm/worklet/WorkletMapField.h>
#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<vtkm::Float64, 3> Center;
vtkm::Float64 Radius;
vtkm::Float64 Omega;
vtkm::Float64 Zeta;
};
}
class OscillatorSource : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(FieldIn<Vec3>, FieldOut<Scalar>);
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<vtkm::Float64, 3>& 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<vtkm::Float64, 3> 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<vtkm::Float64, 3> 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<vtkm::Float64, 3> 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 <typename T>
VTKM_EXEC vtkm::Float64 operator()(const vtkm::Vec<T, 3>& vec) const
{
return (*this)(vtkm::make_Vec(static_cast<vtkm::Float64>(vec[0]),
static_cast<vtkm::Float64>(vec[1]),
static_cast<vtkm::Float64>(vec[2])));
}
private:
vtkm::Vec<internal::Oscillator, MAX_OSCILLATORS> PeriodicOscillators;
vtkm::Vec<internal::Oscillator, MAX_OSCILLATORS> DampedOscillators;
vtkm::Vec<internal::Oscillator, MAX_OSCILLATORS> DecayingOscillators;
vtkm::UInt8 NumberOfPeriodics;
vtkm::UInt8 NumberOfDamped;
vtkm::UInt8 NumberOfDecaying;
vtkm::Float64 Time;
};
}
} // namespace vtkm
#endif // vtk_m_worklet_PointElevation_h