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