diff --git a/examples/moments/CMakeLists.txt b/examples/moments/CMakeLists.txt new file mode 100644 index 000000000..de40a5236 --- /dev/null +++ b/examples/moments/CMakeLists.txt @@ -0,0 +1,35 @@ +##============================================================================= +## +## 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(HelloWorld CXX) + +#Find the VTK-m package +find_package(VTKm REQUIRED QUIET) +vtkm_find_gl(OPTIONAL GL GLUT GLEW) + +#if(TARGET OpenGL::GL AND +# TARGET GLUT::GLUT AND +# TARGET GLEW::GLEW) + add_executable(moments_SERIAL moments.cxx) + target_link_libraries(moments_SERIAL + PRIVATE vtkm_cont ${gl_libs}) +#endif() \ No newline at end of file diff --git a/examples/moments/moments.cxx b/examples/moments/moments.cxx new file mode 100644 index 000000000..f97e46c82 --- /dev/null +++ b/examples/moments/moments.cxx @@ -0,0 +1,151 @@ +// +// Created by ollie on 10/10/18. +// + +#include +#include + +#include +#include + +#include +#include + +namespace detail +{ + +template +vtkm::Float64 eval_functor(vtkm::Float64 x, vtkm::Float64 y); + +template +vtkm::Float64 eval_functor_impl(vtkm::Float64 x, vtkm::Float64 y) +{ + return eval_functor(x, y) * eval_functor(x, y); +} + +template +vtkm::Float64 eval_functor(vtkm::Float64 x, vtkm::Float64 y) +{ + return eval_functor_impl(x, y); +} + +template <> +vtkm::Float64 eval_functor<0>(vtkm::Float64 x, vtkm::Float64 y) +{ + return x; +} + +template <> +vtkm::Float64 eval_functor<1>(vtkm::Float64 x, vtkm::Float64 y) +{ + return y; +} + +template <> +vtkm::Float64 eval_functor<>(vtkm::Float64 x, vtkm::Float64 y) +{ + return 1; +} + +} // detail + +struct Monomial +{ +public: + Monomial(int _p, int _q) + : p(_p) + , q(_q) + { + } + + VTKM_EXEC + vtkm::Float64 eval(vtkm::Float64 x, vtkm::Float64 y) const + { + // return detail::eval_functor(x, y); + return pow(x, p) * pow(y, q); + } + +private: + int p; + int q; +}; + +// FIXME: Why can't we pass the template parameter Radius to WorkletPointNeighborhood? +// TODO: figure out how to make Radius runtime configurable. We need to update +// VTKM's implementation of WorkletPointNeighborhood +template +struct Moments : public vtkm::worklet::WorkletPointNeighborhood<2> +{ +public: + Moments(int _p, int _q) + : monomial(_p, _q) + { + assert(_p >= 0); + assert(_q >= 0); + } + + using ControlSignature = void(CellSetIn, FieldInNeighborhood<>, FieldOut<>); + + using ExecutionSignature = void(_2, _3); + + template + VTKM_EXEC void operator()(const NeighIn& image, vtkm::Float64& moment) const + { + vtkm::Float64 sum = 0; + vtkm::Float64 recp = 1.0 / Radius; + + for (int j = -Radius; j <= Radius; j++) + { + for (int i = -Radius; i <= Radius; i++) + { + if (i * i + j * j <= Radius * Radius) + // TODO: there is several order of magnitude of difference between this and the + // the reference implementation. Scale by dx dy? + sum += monomial.eval(i * recp, j * recp) * image.Get(i, j, 0); + } + } + + moment = sum; + } + +private: + Monomial monomial; +}; + +int main(int argc, char* argv[]) +{ + vtkm::io::reader::VTKDataSetReader reader(argv[1]); + vtkm::cont::DataSet input = reader.ReadDataSet(); + + auto field = input.GetPointField("values"); + vtkm::cont::ArrayHandle pixels; + field.GetData().CopyTo(pixels); + + vtkm::cont::DataSet output; + output.AddCellSet(input.GetCellSet(0)); + output.AddCoordinateSystem(input.GetCoordinateSystem(0)); + + for (int order = 0; order <= 2; ++order) + { + for (int p = order; p >= 0; --p) + { + std::cout << "(" << p << ", " << order - p << ")" << std::endl; + + vtkm::cont::ArrayHandle moments; + + using DispatcherType = vtkm::worklet::DispatcherPointNeighborhood>; + DispatcherType dispatcher(Moments<2>{ p, order - p }); + dispatcher.SetDevice(vtkm::cont::DeviceAdapterTagSerial()); + dispatcher.Invoke(input.GetCellSet(0), pixels, moments); + + std::string fieldName = "moments_"; + fieldName += std::to_string(p) + std::to_string(order - p); + + vtkm::cont::Field momentsField(fieldName, vtkm::cont::Field::Association::POINTS, moments); + output.AddField(momentsField); + } + } + + vtkm::io::writer::VTKDataSetWriter writer("moments_all.vtk"); + writer.WriteDataSet(output); +} \ No newline at end of file