1 //============================================================================ 2 // Copyright (c) Kitware, Inc. 3 // All rights reserved. 4 // See LICENSE.txt for details. 5 // This software is distributed WITHOUT ANY WARRANTY; without even 6 // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 7 // PURPOSE. See the above copyright notice for more information. 8 // 9 // Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS). 10 // Copyright 2014 UT-Battelle, LLC. 11 // Copyright 2014 Los Alamos National Security. 12 // 13 // Under the terms of Contract DE-NA0003525 with NTESS, 14 // the U.S. Government retains certain rights in this software. 15 // 16 // Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National 17 // Laboratory (LANL), the U.S. Government retains certain rights in 18 // this software. 19 //============================================================================ 20 #ifndef vtk_m_worklet_OscillatorSource_h 21 #define vtk_m_worklet_OscillatorSource_h 22 23 #include <vtkm/Math.h> 24 #include <vtkm/worklet/WorkletMapField.h> 25 26 #define MAX_OSCILLATORS 10 27 28 namespace vtkm 29 { 30 namespace worklet 31 { 32 namespace internal 33 { 34 35 struct Oscillator 36 { SetOscillator37 void Set(vtkm::Float64 x, 38 vtkm::Float64 y, 39 vtkm::Float64 z, 40 vtkm::Float64 radius, 41 vtkm::Float64 omega, 42 vtkm::Float64 zeta) 43 { 44 this->Center[0] = x; 45 this->Center[1] = y; 46 this->Center[2] = z; 47 this->Radius = radius; 48 this->Omega = omega; 49 this->Zeta = zeta; 50 } 51 52 vtkm::Vec<vtkm::Float64, 3> Center; 53 vtkm::Float64 Radius; 54 vtkm::Float64 Omega; 55 vtkm::Float64 Zeta; 56 }; 57 } 58 59 class OscillatorSource : public vtkm::worklet::WorkletMapField 60 { 61 public: 62 typedef void ControlSignature(FieldIn<Vec3>, FieldOut<Scalar>); 63 typedef _2 ExecutionSignature(_1); 64 65 VTKM_CONT OscillatorSource()66 OscillatorSource() 67 : NumberOfPeriodics(0) 68 , NumberOfDamped(0) 69 , NumberOfDecaying(0) 70 { 71 } 72 73 VTKM_CONT AddPeriodic(vtkm::Float64 x,vtkm::Float64 y,vtkm::Float64 z,vtkm::Float64 radius,vtkm::Float64 omega,vtkm::Float64 zeta)74 void AddPeriodic(vtkm::Float64 x, 75 vtkm::Float64 y, 76 vtkm::Float64 z, 77 vtkm::Float64 radius, 78 vtkm::Float64 omega, 79 vtkm::Float64 zeta) 80 { 81 if (this->NumberOfPeriodics < MAX_OSCILLATORS) 82 { 83 this->PeriodicOscillators[this->NumberOfPeriodics].Set(x, y, z, radius, omega, zeta); 84 this->NumberOfPeriodics++; 85 } 86 } 87 88 VTKM_CONT AddDamped(vtkm::Float64 x,vtkm::Float64 y,vtkm::Float64 z,vtkm::Float64 radius,vtkm::Float64 omega,vtkm::Float64 zeta)89 void AddDamped(vtkm::Float64 x, 90 vtkm::Float64 y, 91 vtkm::Float64 z, 92 vtkm::Float64 radius, 93 vtkm::Float64 omega, 94 vtkm::Float64 zeta) 95 { 96 if (this->NumberOfDamped < MAX_OSCILLATORS) 97 { 98 this->DampedOscillators[this->NumberOfDamped * 6].Set(x, y, z, radius, omega, zeta); 99 this->NumberOfDamped++; 100 } 101 } 102 103 VTKM_CONT AddDecaying(vtkm::Float64 x,vtkm::Float64 y,vtkm::Float64 z,vtkm::Float64 radius,vtkm::Float64 omega,vtkm::Float64 zeta)104 void AddDecaying(vtkm::Float64 x, 105 vtkm::Float64 y, 106 vtkm::Float64 z, 107 vtkm::Float64 radius, 108 vtkm::Float64 omega, 109 vtkm::Float64 zeta) 110 { 111 if (this->NumberOfDecaying < MAX_OSCILLATORS) 112 { 113 this->DecayingOscillators[this->NumberOfDecaying * 6].Set(x, y, z, radius, omega, zeta); 114 this->NumberOfDecaying++; 115 } 116 } 117 118 VTKM_CONT SetTime(vtkm::Float64 time)119 void SetTime(vtkm::Float64 time) { this->Time = time; } 120 121 VTKM_EXEC operator()122 vtkm::Float64 operator()(const vtkm::Vec<vtkm::Float64, 3>& vec) const 123 { 124 vtkm::UInt8 oIdx; 125 vtkm::Float64 t0, t, result = 0; 126 const internal::Oscillator* oscillator; 127 128 t0 = 0.0; 129 t = this->Time * 2 * 3.14159265358979323846; 130 131 // Compute damped 132 for (oIdx = 0; oIdx < this->NumberOfDamped; oIdx++) 133 { 134 oscillator = &this->DampedOscillators[oIdx]; 135 136 vtkm::Vec<vtkm::Float64, 3> delta = oscillator->Center - vec; 137 vtkm::Float64 dist2 = dot(delta, delta); 138 vtkm::Float64 dist_damp = vtkm::Exp(-dist2 / (2 * oscillator->Radius * oscillator->Radius)); 139 vtkm::Float64 phi = vtkm::ACos(oscillator->Zeta); 140 vtkm::Float64 val = 1. - 141 vtkm::Exp(-oscillator->Zeta * oscillator->Omega * t0) * 142 (vtkm::Sin(vtkm::Sqrt(1 - oscillator->Zeta * oscillator->Zeta) * oscillator->Omega * t + 143 phi) / 144 vtkm::Sin(phi)); 145 result += val * dist_damp; 146 } 147 148 // Compute decaying 149 for (oIdx = 0; oIdx < this->NumberOfDecaying; oIdx++) 150 { 151 oscillator = &this->DecayingOscillators[oIdx]; 152 t = t0 + 1 / oscillator->Omega; 153 vtkm::Vec<vtkm::Float64, 3> delta = oscillator->Center - vec; 154 vtkm::Float64 dist2 = dot(delta, delta); 155 vtkm::Float64 dist_damp = vtkm::Exp(-dist2 / (2 * oscillator->Radius * oscillator->Radius)); 156 vtkm::Float64 val = vtkm::Sin(t / oscillator->Omega) / (oscillator->Omega * t); 157 result += val * dist_damp; 158 } 159 160 // Compute periodic 161 for (oIdx = 0; oIdx < this->NumberOfPeriodics; oIdx++) 162 { 163 oscillator = &this->PeriodicOscillators[oIdx]; 164 t = t0 + 1 / oscillator->Omega; 165 vtkm::Vec<vtkm::Float64, 3> delta = oscillator->Center - vec; 166 vtkm::Float64 dist2 = dot(delta, delta); 167 vtkm::Float64 dist_damp = vtkm::Exp(-dist2 / (2 * oscillator->Radius * oscillator->Radius)); 168 vtkm::Float64 val = vtkm::Sin(t / oscillator->Omega); 169 result += val * dist_damp; 170 } 171 172 // We are done... 173 return result; 174 } 175 176 template <typename T> operator()177 VTKM_EXEC vtkm::Float64 operator()(const vtkm::Vec<T, 3>& vec) const 178 { 179 return (*this)(vtkm::make_Vec(static_cast<vtkm::Float64>(vec[0]), 180 static_cast<vtkm::Float64>(vec[1]), 181 static_cast<vtkm::Float64>(vec[2]))); 182 } 183 184 private: 185 vtkm::Vec<internal::Oscillator, MAX_OSCILLATORS> PeriodicOscillators; 186 vtkm::Vec<internal::Oscillator, MAX_OSCILLATORS> DampedOscillators; 187 vtkm::Vec<internal::Oscillator, MAX_OSCILLATORS> DecayingOscillators; 188 vtkm::UInt8 NumberOfPeriodics; 189 vtkm::UInt8 NumberOfDamped; 190 vtkm::UInt8 NumberOfDecaying; 191 vtkm::Float64 Time; 192 }; 193 } 194 195 } // namespace vtkm 196 197 #endif // vtk_m_worklet_PointElevation_h 198