UnitTestCellGradient.cxx 10.9 KB
Newer Older
1 2 3 4 5 6 7 8
//============================================================================
//  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.
//
Kenneth Moreland's avatar
Kenneth Moreland committed
9
//  Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
10 11 12
//  Copyright 2014 UT-Battelle, LLC.
//  Copyright 2014 Los Alamos National Security.
//
Kenneth Moreland's avatar
Kenneth Moreland committed
13
//  Under the terms of Contract DE-NA0003525 with NTESS,
14 15 16 17 18 19 20
//  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.
//============================================================================

21
#include <vtkm/worklet/Gradient.h>
22 23

#include <vtkm/cont/testing/MakeTestDataSet.h>
24
#include <vtkm/cont/testing/Testing.h>
25

26 27
namespace
{
28

29
template <typename DeviceAdapter>
30 31 32 33 34 35 36
void TestCellGradientUniform2D()
{
  std::cout << "Testing CellGradient Worklet on 2D structured data" << std::endl;

  vtkm::cont::testing::MakeTestDataSet testDataSet;
  vtkm::cont::DataSet dataSet = testDataSet.Make2DUniformDataSet0();

37
  vtkm::cont::ArrayHandle<vtkm::Float32> input;
38
  vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 3>> result;
39

40 41 42 43 44
  dataSet.GetField("pointvar").GetData().CopyTo(input);

  vtkm::worklet::CellGradient gradient;
  result =
    gradient.Run(dataSet.GetCellSet(), dataSet.GetCoordinateSystem(), input, DeviceAdapter());
45

46
  vtkm::Vec<vtkm::Float32, 3> expected[2] = { { 10, 30, 0 }, { 10, 30, 0 } };
47 48
  for (int i = 0; i < 2; ++i)
  {
49 50
    VTKM_TEST_ASSERT(test_equal(result.GetPortalConstControl().Get(i), expected[i]),
                     "Wrong result for CellGradient worklet on 2D uniform data");
51 52 53
  }
}

54
template <typename DeviceAdapter>
55 56
void TestCellGradientUniform3D()
{
luz.paz's avatar
luz.paz committed
57
  std::cout << "Testing CellGradient Worklet on 3D structured data" << std::endl;
58 59 60 61

  vtkm::cont::testing::MakeTestDataSet testDataSet;
  vtkm::cont::DataSet dataSet = testDataSet.Make3DUniformDataSet0();

62 63 64 65
  vtkm::cont::ArrayHandle<vtkm::Float32> input;
  vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 3>> result;

  dataSet.GetField("pointvar").GetData().CopyTo(input);
66

67 68 69
  vtkm::worklet::CellGradient gradient;
  result =
    gradient.Run(dataSet.GetCellSet(), dataSet.GetCoordinateSystem(), input, DeviceAdapter());
70

71
  vtkm::Vec<vtkm::Float32, 3> expected[4] = {
72 73 74 75
    { 10.025f, 30.075f, 60.125f },
    { 10.025f, 30.075f, 60.125f },
    { 10.025f, 30.075f, 60.175f },
    { 10.025f, 30.075f, 60.175f },
76
  };
77 78
  for (int i = 0; i < 4; ++i)
  {
79 80
    VTKM_TEST_ASSERT(test_equal(result.GetPortalConstControl().Get(i), expected[i]),
                     "Wrong result for CellGradient worklet on 3D uniform data");
81
  }
82 83
}

84
template <typename DeviceAdapter>
85 86
void TestCellGradientUniform3DWithVectorField()
{
87
  std::cout
luz.paz's avatar
luz.paz committed
88
    << "Testing CellGradient and QCriterion Worklet with a vector field on 3D structured data"
89
    << std::endl;
90 91 92 93 94
  vtkm::cont::testing::MakeTestDataSet testDataSet;
  vtkm::cont::DataSet dataSet = testDataSet.Make3DUniformDataSet0();

  //Verify that we can compute the gradient of a 3 component vector
  const int nVerts = 18;
95 96 97 98
  vtkm::Float64 vars[nVerts] = { 10.1,  20.1,  30.1,  40.1,  50.2,  60.2,  70.2,  80.2,  90.3,
                                 100.3, 110.3, 120.3, 130.4, 140.4, 150.4, 160.4, 170.5, 180.5 };
  std::vector<vtkm::Vec<vtkm::Float64, 3>> vec(18);
  for (std::size_t i = 0; i < vec.size(); ++i)
99
  {
100
    vec[i] = vtkm::make_Vec(vars[i], vars[i], vars[i]);
101
  }
102
  vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float64, 3>> input = vtkm::cont::make_ArrayHandle(vec);
103 104

  //we need to add Vec3 array to the dataset
105
  vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Vec<vtkm::Float64, 3>, 3>> result;
106 107 108 109 110 111

  vtkm::worklet::GradientOutputFields<vtkm::Vec<vtkm::Float64, 3>> extraOutput;
  extraOutput.SetComputeDivergence(false);
  extraOutput.SetComputeVorticity(false);
  extraOutput.SetComputeQCriterion(true);

112
  vtkm::worklet::CellGradient gradient;
113 114 115 116 117 118 119 120 121 122 123
  result = gradient.Run(
    dataSet.GetCellSet(), dataSet.GetCoordinateSystem(), input, extraOutput, DeviceAdapter());

  VTKM_TEST_ASSERT((extraOutput.Gradient.GetNumberOfValues() == 4),
                   "Gradient field should be generated");
  VTKM_TEST_ASSERT((extraOutput.Divergence.GetNumberOfValues() == 0),
                   "Divergence field shouldn't be generated");
  VTKM_TEST_ASSERT((extraOutput.Vorticity.GetNumberOfValues() == 0),
                   "Vorticity field shouldn't be generated");
  VTKM_TEST_ASSERT((extraOutput.QCriterion.GetNumberOfValues() == 4),
                   "QCriterion field should be generated");
124 125 126 127 128 129 130

  vtkm::Vec<vtkm::Vec<vtkm::Float64, 3>, 3> expected[4] = {
    { { 10.025, 10.025, 10.025 }, { 30.075, 30.075, 30.075 }, { 60.125, 60.125, 60.125 } },
    { { 10.025, 10.025, 10.025 }, { 30.075, 30.075, 30.075 }, { 60.125, 60.125, 60.125 } },
    { { 10.025, 10.025, 10.025 }, { 30.075, 30.075, 30.075 }, { 60.175, 60.175, 60.175 } },
    { { 10.025, 10.025, 10.025 }, { 30.075, 30.075, 30.075 }, { 60.175, 60.175, 60.175 } }
  };
131 132
  for (int i = 0; i < 4; ++i)
  {
133 134 135 136 137 138 139 140 141
    vtkm::Vec<vtkm::Vec<vtkm::Float64, 3>, 3> e = expected[i];
    vtkm::Vec<vtkm::Vec<vtkm::Float64, 3>, 3> r = result.GetPortalConstControl().Get(i);

    VTKM_TEST_ASSERT(test_equal(e[0], r[0]),
                     "Wrong result for vec field CellGradient worklet on 3D uniform data");
    VTKM_TEST_ASSERT(test_equal(e[1], r[1]),
                     "Wrong result for vec field CellGradient worklet on 3D uniform data");
    VTKM_TEST_ASSERT(test_equal(e[2], r[2]),
                     "Wrong result for vec field CellGradient worklet on 3D uniform data");
142

143 144
    const vtkm::Vec<vtkm::Float64, 3> v(e[1][2] - e[2][1], e[2][0] - e[0][2], e[0][1] - e[1][0]);
    const vtkm::Vec<vtkm::Float64, 3> s(e[1][2] + e[2][1], e[2][0] + e[0][2], e[0][1] + e[1][0]);
145 146 147 148
    const vtkm::Vec<vtkm::Float64, 3> d(e[0][0], e[1][1], e[2][2]);

    //compute QCriterion
    vtkm::Float64 qcriterion =
149
      ((vtkm::Dot(v, v) / 2.0f) - (vtkm::Dot(d, d) + (vtkm::Dot(s, s) / 2.0f))) / 2.0f;
150 151 152

    vtkm::Float64 q = extraOutput.QCriterion.GetPortalConstControl().Get(i);

153 154 155
    std::cout << "qcriterion expected: " << qcriterion << std::endl;
    std::cout << "qcriterion actual: " << q << std::endl;

156 157 158 159 160 161 162 163 164
    VTKM_TEST_ASSERT(
      test_equal(qcriterion, q),
      "Wrong result for QCriterion field of CellGradient worklet on 3D uniform data");
  }
}

template <typename DeviceAdapter>
void TestCellGradientUniform3DWithVectorField2()
{
luz.paz's avatar
luz.paz committed
165
  std::cout << "Testing CellGradient Worklet with a vector field on 3D structured data" << std::endl
166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218
            << "Disabling Gradient computation and enabling Divergence, and Vorticity" << std::endl;
  vtkm::cont::testing::MakeTestDataSet testDataSet;
  vtkm::cont::DataSet dataSet = testDataSet.Make3DUniformDataSet0();

  //Verify that we can compute the gradient of a 3 component vector
  const int nVerts = 18;
  vtkm::Float64 vars[nVerts] = { 10.1,  20.1,  30.1,  40.1,  50.2,  60.2,  70.2,  80.2,  90.3,
                                 100.3, 110.3, 120.3, 130.4, 140.4, 150.4, 160.4, 170.5, 180.5 };
  std::vector<vtkm::Vec<vtkm::Float64, 3>> vec(18);
  for (std::size_t i = 0; i < vec.size(); ++i)
  {
    vec[i] = vtkm::make_Vec(vars[i], vars[i], vars[i]);
  }
  vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float64, 3>> input = vtkm::cont::make_ArrayHandle(vec);

  vtkm::worklet::GradientOutputFields<vtkm::Vec<vtkm::Float64, 3>> extraOutput;
  extraOutput.SetComputeGradient(false);
  extraOutput.SetComputeDivergence(true);
  extraOutput.SetComputeVorticity(true);
  extraOutput.SetComputeQCriterion(false);

  vtkm::worklet::CellGradient gradient;
  auto result = gradient.Run(
    dataSet.GetCellSet(), dataSet.GetCoordinateSystem(), input, extraOutput, DeviceAdapter());

  //Verify that the result is 0 size
  VTKM_TEST_ASSERT((result.GetNumberOfValues() == 0), "Gradient field shouldn't be generated");
  //Verify that the extra arrays are the correct size
  VTKM_TEST_ASSERT((extraOutput.Gradient.GetNumberOfValues() == 0),
                   "Gradient field shouldn't be generated");
  VTKM_TEST_ASSERT((extraOutput.Divergence.GetNumberOfValues() == 4),
                   "Divergence field should be generated");
  VTKM_TEST_ASSERT((extraOutput.Vorticity.GetNumberOfValues() == 4),
                   "Vorticity field should be generated");
  VTKM_TEST_ASSERT((extraOutput.QCriterion.GetNumberOfValues() == 0),
                   "QCriterion field shouldn't be generated");

  //Verify the contents of the other arrays
  vtkm::Vec<vtkm::Vec<vtkm::Float64, 3>, 3> expected_gradients[4] = {
    { { 10.025, 10.025, 10.025 }, { 30.075, 30.075, 30.075 }, { 60.125, 60.125, 60.125 } },
    { { 10.025, 10.025, 10.025 }, { 30.075, 30.075, 30.075 }, { 60.125, 60.125, 60.125 } },
    { { 10.025, 10.025, 10.025 }, { 30.075, 30.075, 30.075 }, { 60.175, 60.175, 60.175 } },
    { { 10.025, 10.025, 10.025 }, { 30.075, 30.075, 30.075 }, { 60.175, 60.175, 60.175 } }
  };

  for (int i = 0; i < 4; ++i)
  {
    vtkm::Vec<vtkm::Vec<vtkm::Float64, 3>, 3> eg = expected_gradients[i];

    vtkm::Float64 d = extraOutput.Divergence.GetPortalConstControl().Get(i);
    VTKM_TEST_ASSERT(test_equal((eg[0][0] + eg[1][1] + eg[2][2]), d),
                     "Wrong result for Divergence on 3D uniform data");

219
    vtkm::Vec<vtkm::Float64, 3> ev(eg[1][2] - eg[2][1], eg[2][0] - eg[0][2], eg[0][1] - eg[1][0]);
220 221
    vtkm::Vec<vtkm::Float64, 3> v = extraOutput.Vorticity.GetPortalConstControl().Get(i);
    VTKM_TEST_ASSERT(test_equal(ev, v), "Wrong result for Vorticity on 3D uniform data");
222
  }
223 224
}

225
template <typename DeviceAdapter>
226 227 228 229 230 231 232
void TestCellGradientExplicit()
{
  std::cout << "Testing CellGradient Worklet on Explicit data" << std::endl;

  vtkm::cont::testing::MakeTestDataSet testDataSet;
  vtkm::cont::DataSet dataSet = testDataSet.Make3DExplicitDataSet0();

233
  vtkm::cont::ArrayHandle<vtkm::Float32> input;
234
  vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 3>> result;
235
  dataSet.GetField("pointvar").GetData().CopyTo(input);
236

237 238 239
  vtkm::worklet::CellGradient gradient;
  result =
    gradient.Run(dataSet.GetCellSet(), dataSet.GetCoordinateSystem(), input, DeviceAdapter());
240

241
  vtkm::Vec<vtkm::Float32, 3> expected[2] = { { 10.f, 10.1f, 0.0f }, { 10.f, 10.1f, -0.0f } };
242 243
  for (int i = 0; i < 2; ++i)
  {
244 245
    VTKM_TEST_ASSERT(test_equal(result.GetPortalConstControl().Get(i), expected[i]),
                     "Wrong result for CellGradient worklet on 3D explicit data");
246 247 248 249 250
  }
}

void TestCellGradient()
{
251 252 253 254
  using DeviceAdapter = VTKM_DEFAULT_DEVICE_ADAPTER_TAG;
  TestCellGradientUniform2D<DeviceAdapter>();
  TestCellGradientUniform3D<DeviceAdapter>();
  TestCellGradientUniform3DWithVectorField<DeviceAdapter>();
255
  TestCellGradientUniform3DWithVectorField2<DeviceAdapter>();
256
  TestCellGradientExplicit<DeviceAdapter>();
257 258 259
}
}

260
int UnitTestCellGradient(int, char* [])
261 262 263
{
  return vtkm::cont::testing::Testing::Run(TestCellGradient);
}