Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-19 18:45:43 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Change GaussianSplatter to KernelSplatter to support other kernels

Browse Source 
Template the splatter algorithm over Kernel type and use abstract kernel
interface to fetch the kernel value.

Add Gaussian, Spline3rdOrder kernel classes templated over dimension.
Other kernels can/will be added in future.

Kernel classes are defined such that the integral of the volume is unity
as is the convention in (for example) SPH simulations.
This commit is contained in:
John Biddiscombe 2015-09-14 09:27:32 +02:00
parent 378cb17e8a
commit 29001e377f
4 changed files with 547 additions and 66 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

136
vtkm/worklet/GaussianSplatter.h → vtkm/worklet/KernelSplatter.h
View File

@ -17,8 +17,8 @@
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_GaussianSplatter_h
#define vtk_m_worklet_GaussianSplatter_h
#ifndef vtk_m_worklet_KernelSplatter_h
#define vtk_m_worklet_KernelSplatter_h
#include <vtkm/Math.h>
@ -38,6 +38,10 @@
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/kernels/KernelBase.h>
#include <vtkm/worklet/kernels/Gaussian.h>
#include <vtkm/worklet/kernels/Spline3rdOrder.h>
#define __VTKM_GAUSSIAN_SPLATTER_BENCHMARK
//----------------------------------------------------------------------------
@ -57,7 +61,7 @@
#endif
//----------------------------------------------------------------------------
// Gaussian splatter worklet/filter
// Kernel splatter worklet/filter
//----------------------------------------------------------------------------
namespace vtkm { namespace worklet
{
@ -123,8 +127,8 @@ void OutputArrayDebug(const vtkm::cont::ArrayHandlePermutation<I, vtkm::cont::Ar
#endif
} // namespace debug
template<typename DeviceAdapter>
struct GaussianSplatterFilterUniformGrid
template<typename Kernel, typename DeviceAdapter>
struct KernelSplatterFilterUniformGrid
{
typedef vtkm::cont::ArrayHandle<vtkm::Float64> DoubleHandleType;
typedef vtkm::cont::ArrayHandle<vtkm::Float32> FloatHandleType;
@ -170,9 +174,10 @@ struct GaussianSplatterFilterUniformGrid
class GetFootprint : public vtkm::worklet::WorkletMapField
{
private:
vtkm::Vec<vtkm::Float64, 3> Origin;
vtkm::Vec<vtkm::Float64, 3> Spacing;
vtkm::Vec<vtkm::Float64, 3> origin_;
vtkm::Vec<vtkm::Float64, 3> spacing_;
vtkm::Id3 VolumeDimensions;
Kernel kernel_;
public:
typedef void ControlSignature(FieldIn<>, FieldIn<>, FieldIn<>, FieldIn<>,
@ -183,30 +188,32 @@ struct GaussianSplatterFilterUniformGrid
GetFootprint(
const vtkm::Vec<vtkm::Float64, 3> &o,
const vtkm::Vec<vtkm::Float64, 3> &s,
const vtkm::Id3 &dim)
: Origin(o), Spacing(s),
VolumeDimensions(dim) { }
const vtkm::Id3 &dim,
const Kernel &kernel)
: origin_(o), spacing_(s),
VolumeDimensions(dim), kernel_(kernel) { }
template<typename T, typename T2>
VTKM_EXEC_CONT_EXPORT
void operator()(
const T &xValue,
const T &yValue,
const T &zValue,
const T2 &rValue,
const T &x,
const T &y,
const T &z,
const T2 &h,
vtkm::Vec<vtkm::Float64, 3> &splatPoint,
vtkm::Id3 &minFootprint,
vtkm::Id3 &maxFootprint,
vtkm::Id &footprintSize) const
{
PointType splat, min, max;
vtkm::Vec<vtkm::Float64, 3> sample = vtkm::make_Vec(xValue, yValue, zValue);
vtkm::Vec<vtkm::Float64, 3> sample = vtkm::make_Vec(x, y, z);
vtkm::Id size = 1;
double cutoff = kernel_.maxDistance(h);
for (int i = 0; i < 3; i++)
{
splat[i] = (sample[i] - this->Origin[i]) / this->Spacing[i];
min[i] = static_cast<vtkm::Id>(ceil(static_cast<double>(splat[i]) - rValue));
max[i] = static_cast<vtkm::Id>(floor(static_cast<double>(splat[i]) + rValue));
splat[i] = (sample[i] - this->origin_[i]) / this->spacing_[i];
min[i] = static_cast<vtkm::Id>(ceil(static_cast<double>(splat[i]) - cutoff));
max[i] = static_cast<vtkm::Id>(floor(static_cast<double>(splat[i]) + cutoff));
if (min[i] < 0) {
min[i] = 0;
}
@ -246,35 +253,34 @@ struct GaussianSplatterFilterUniformGrid
};
//-----------------------------------------------------------------------
// Compute the Gaussian splat value of the input neighbour point.
// Compute the splat value of the input neighbour point.
// The voxel Id of this point within the volume is also determined.
// @TODO : fix kernel to behave with correct radious, scale exponent etc
//-----------------------------------------------------------------------
class GetSplatValue : public vtkm::worklet::WorkletMapField
{
private:
vtkm::Vec<vtkm::Float64, 3> Spacing;
vtkm::Vec<vtkm::Float64, 3> Origin;
vtkm::Vec<vtkm::Float64, 3> spacing_;
vtkm::Vec<vtkm::Float64, 3> origin_;
vtkm::Id3 VolumeDim;
vtkm::Float64 Radius2;
vtkm::Float64 ExponentFactor;
vtkm::Float64 ScalingFactor;
Kernel kernel;
public:
typedef void ControlSignature(FieldIn<>, FieldIn<>, FieldIn<>, FieldIn<>,
typedef void ControlSignature(
FieldIn<>, FieldIn<>, FieldIn<>, FieldIn<>, FieldIn<>,
FieldIn<>, FieldIn<>, FieldOut<>,
FieldOut<>);
typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6, _7, _8);
typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6, _7, _8, _9);
VTKM_CONT_EXPORT
GetSplatValue(
const vtkm::Vec<vtkm::Float64, 3> &orig,
const vtkm::Vec<vtkm::Float64, 3> &s,
const vtkm::Id3 &dim,
const vtkm::Float64 &ef,
const vtkm::Float64 &sf)
: Spacing(s), Origin(orig), VolumeDim(dim),
ExponentFactor(ef), ScalingFactor(sf) {}
const Kernel &k)
: spacing_(s), origin_(orig), VolumeDim(dim), kernel(k) {}
template<typename T, typename T2, typename P>
VTKM_EXEC_CONT_EXPORT
@ -282,7 +288,8 @@ struct GaussianSplatterFilterUniformGrid
const vtkm::Vec<P, 3> &splatPoint,
const T &minBound,
const T &maxBound,
const T2 &radius,
const T2 &kernel_H,
const T2 &scale,
const vtkm::Id splatPointId,
const vtkm::Id localNeighborId,
vtkm::Id &neighborVoxelId,
@ -295,23 +302,18 @@ struct GaussianSplatterFilterUniformGrid
vtkm::Id remainder = localNeighborId % divisor;
vtkm::Id j = remainder / xRange;
vtkm::Id k = remainder % xRange;
vtkm::Float64 radius2 = radius*radius;
vtkm::Id3 voxel = minBound + vtkm::make_Vec(k, j, i);
// note the order of k,j,i
vtkm::Id3 voxel = minBound + vtkm::make_Vec(k, j, i);
PointType dist = vtkm::make_Vec(
(splatPoint[0] - voxel[0])*Spacing[0],
(splatPoint[1] - voxel[1])*Spacing[0],
(splatPoint[2] - voxel[2])*Spacing[0]
(splatPoint[0] - voxel[0])*spacing_[0],
(splatPoint[1] - voxel[1])*spacing_[0],
(splatPoint[2] - voxel[2])*spacing_[0]
);
vtkm::Float64 dist2 = vtkm::dot(dist,dist);
// Compute Gaussian splat value
if (dist2 <= radius2) {
splatValue = this->ScalingFactor *
vtkm::Exp((this->ExponentFactor*(dist2) / radius2));
}
else {
splatValue = 0.0;
}
// Compute splat value using the kernel distance_squared function
splatValue = scale*kernel.w2(kernel_H, dist2);
//
neighborVoxelId = (voxel[2] * VolumeDim[0] * VolumeDim[1]) +
(voxel[1] * VolumeDim[0]) + voxel[0];
if (neighborVoxelId<0) neighborVoxelId = -1;
@ -333,7 +335,6 @@ struct GaussianSplatterFilterUniformGrid
VTKM_CONT_EXPORT
UpdateVoxelSplats() {}
VTKM_EXEC_CONT_EXPORT
void operator()(const vtkm::Id &voxelIndex,
const vtkm::Float64 &splatValue,
@ -346,26 +347,31 @@ struct GaussianSplatterFilterUniformGrid
//-----------------------------------------------------------------------
// Construct a splatter filter/object
//
// @TODO, get the origin and spacing from the dataset coordinates
// @TODO, get the origin_ and spacing_ from the dataset coordinates
// instead of requiring them to be passed as parameters.
//-----------------------------------------------------------------------
GaussianSplatterFilterUniformGrid(const vtkm::Id3 &dims,
KernelSplatterFilterUniformGrid(const vtkm::Id3 &dims,
vtkm::Vec<vtkm::FloatDefault, 3> origin,
vtkm::Vec<vtkm::FloatDefault, 3> spacing,
const vtkm::cont::DataSet &dataSet) :
CDims(dims),
Origin(origin),
Spacing(spacing),
DataSet(dataSet)
const vtkm::cont::DataSet &dataset,
const Kernel &kernel) :
dims_(dims),
origin_(origin),
spacing_(spacing),
dataset_(dataset),
kernel_(kernel)
{
}
//-----------------------------------------------------------------------
// class variables for the splat filter
//-----------------------------------------------------------------------
vtkm::Id3 CDims;
FloatVec Origin, Spacing;
vtkm::cont::DataSet DataSet;
vtkm::Id3 dims_;
FloatVec origin_;
FloatVec spacing_;
vtkm::cont::DataSet dataset_;
// The kernel used for this filter
Kernel kernel_;
//-----------------------------------------------------------------------
// Run the filter, given the input params
@ -376,15 +382,12 @@ struct GaussianSplatterFilterUniformGrid
const vtkm::cont::ArrayHandle<vtkm::Float64, StorageT> yValues,
const vtkm::cont::ArrayHandle<vtkm::Float64, StorageT> zValues,
const vtkm::cont::ArrayHandle<vtkm::Float32, StorageT> rValues,
const vtkm::cont::ArrayHandle<vtkm::Float32, StorageT> sValues,
FloatHandleType scalarSplatOutput)
{
// Define the constants for the algorithm
const vtkm::Float64 exponentFactor = -5.0;
const vtkm::Float64 scaleFactor = 1.0;
// Number of grid points in the volume bounding box
vtkm::Id3 pointDimensions = vtkm::make_Vec(CDims[0]+1, CDims[1]+1, CDims[2]+1);
const vtkm::Id numVolumePoints = (CDims[0] + 1) * (CDims[1] + 1) * (CDims[2] + 1);
vtkm::Id3 pointDimensions = vtkm::make_Vec(dims_[0]+1, dims_[1]+1, dims_[2]+1);
const vtkm::Id numVolumePoints = (dims_[0] + 1) * (dims_[1] + 1) * (dims_[2] + 1);
//---------------------------------------------------------------
// Get the splat footprint/neighborhood of each sample point, as
@ -396,7 +399,7 @@ struct GaussianSplatterFilterUniformGrid
IdHandleType numNeighbors;
IdHandleType localNeighborIds;
GetFootprint footprint_worklet(Origin, Spacing, pointDimensions);
GetFootprint footprint_worklet(origin_, spacing_, pointDimensions, kernel_);
vtkm::worklet::DispatcherMapField<GetFootprint> footprintDispatcher(footprint_worklet);
START_TIMER_BLOCK(GetFootprint)
@ -489,6 +492,7 @@ struct GaussianSplatterFilterUniformGrid
VecPermType ptFootprintMins(neighbor2SplatId, footprintMin);
VecPermType ptFootprintMaxs(neighbor2SplatId, footprintMax);
FloatPermType radii(neighbor2SplatId, rValues);
FloatPermType scale(neighbor2SplatId, sValues);
debug::OutputArrayDebug(radii, "radii");
debug::OutputArrayDebug(ptSplatPoints, "ptSplatPoints");
@ -503,11 +507,10 @@ struct GaussianSplatterFilterUniformGrid
FloatHandleType splatValues;
GetSplatValue splatterDispatcher_worklet(
Origin,
Spacing,
origin_,
spacing_,
pointDimensions,
exponentFactor,
scaleFactor);
kernel_);
vtkm::worklet::DispatcherMapField<GetSplatValue> splatterDispatcher(splatterDispatcher_worklet);
START_TIMER_BLOCK(GetSplatValue)
@ -516,6 +519,7 @@ struct GaussianSplatterFilterUniformGrid
ptFootprintMins,
ptFootprintMaxs,
radii,
scale,
neighbor2SplatId,
localNeighborIds, neighborVoxelIds, splatValues);
END_TIMER_BLOCK(GetSplatValue)
@ -589,10 +593,10 @@ struct GaussianSplatterFilterUniformGrid
voxelSplatSums.ReleaseResources();
}
}; //struct GaussianSplatter
}; //struct KernelSplatter
}
} //namespace vtkm::worklet
#endif //vtk_m_worklet_GaussianSplatter_h
#endif //vtk_m_worklet_KernelSplatter_h

168
vtkm/worklet/kernels/Gaussian.h Normal file
View File

@ -0,0 +1,168 @@
#ifndef VTKM_KERNEL_GAUSSIAN_H
#define VTKM_KERNEL_GAUSSIAN_H
#include "KernelBase.h"
//
// Gaussian kernel.
// Compact support is achived by truncating the kernel beyond the cutoff radius
// This implementation uses a factor of 5 between smoothing length and cutoff
//
namespace vtkm { namespace worklet {
namespace kernels {
template <int Dimensions>
struct Gaussian : public KernelBase< Gaussian<Dimensions> >
{
//---------------------------------------------------------------------
// Constructor
// Calculate coefficients used repeatedly when evaluating the kernel
// value or gradient
Gaussian(double smoothingLength)
: KernelBase< Gaussian<Dimensions> >(smoothingLength)
{
Hinverse_ = 1.0/smoothingLength;
Hinverse2_ = Hinverse_*Hinverse_;
maxRadius_ = 5.0*smoothingLength;
maxRadius2_ = maxRadius_*maxRadius_;
//
norm_ = 1.0 / pow(M_PI, static_cast<double>(Dimensions) / 2.0);
scale_W_ = norm_ * power<Dimensions> (Hinverse_);
scale_GradW_ = - 2.0 * power<Dimensions+1>(Hinverse_) / norm_;
}
//---------------------------------------------------------------------
// return the multiplier between smoothing length and max cutoff distance
constexpr double getDilationFactor() const { return 5.0; }
//---------------------------------------------------------------------
// compute w(h) for the given distance
inline double w(double distance) const
{
if (distance<maxDistance()) {
// compute r/h
double normedDist = distance * Hinverse_;
// compute w(h)
return scale_W_ * exp(-normedDist * normedDist);
}
return 0.0;
}
//---------------------------------------------------------------------
// compute w(h) for the given squared distance
inline double w2(double distance2) const
{
if (distance2<maxSquaredDistance()) {
// compute (r/h)^2
double normedDist = distance2 * Hinverse2_;
// compute w(h)
return scale_W_ * exp(-normedDist);
}
return 0.0;
}
//---------------------------------------------------------------------
// compute w(h) for a variable h kernel
inline double w(double h, double distance) const
{
if (distance<maxDistance(h)) {
double Hinverse = 1.0/h;
double scale_W = norm_ * power<Dimensions>(Hinverse);
double Q = distance * Hinverse;
return scale_W * exp(-Q*Q);
}
return 0;
}
//---------------------------------------------------------------------
// compute w(h) for a variable h kernel using distance squared
inline double w2(double h, double distance2) const
{
if (distance2<maxSquaredDistance(h)) {
double Hinverse = 1.0/h;
double scale_W = norm_ * power<Dimensions>(Hinverse);
double Q = distance2 * Hinverse * Hinverse;
return scale_W * exp(-Q);
}
return 0;
}
//---------------------------------------------------------------------
// Calculates the kernel derivative for a distance {x,y,z} vector
// from the centre
inline vector_type gradW(double distance, const vector_type& pos) const
{
double Q = distance * Hinverse_;
if (Q != 0.0)
{
return scale_GradW_ * exp(-Q * Q) * pos;
}
else {
return vector_type(0.0);
}
}
//---------------------------------------------------------------------
// Calculates the kernel derivative for a distance {x,y,z} vector
// from the centre using a variable h
inline vector_type gradW(double h, double distance, const vector_type& pos) const
{
double Hinverse = 1.0/h;
double scale_GradW = - 2.0 * power<Dimensions+1>(Hinverse)
/ pow(M_PI, static_cast<double>(Dimensions) / 2.0);
double Q = distance * Hinverse;
//!!! check this due to the fitting offset
if (distance != 0.0)
{
return scale_GradW * exp(-Q * Q) * pos;
}
else {
return vector_type(0.0);
}
}
//---------------------------------------------------------------------
// return the maximum distance at which this kernel is non zero
inline double maxDistance() const
{
return maxRadius_;
}
//---------------------------------------------------------------------
// return the maximum distance at which this variable h kernel is non zero
inline double maxDistance(double h) const
{
return getDilationFactor()*h;
}
//---------------------------------------------------------------------
// return the maximum distance at which this kernel is non zero
inline double maxSquaredDistance() const
{
return maxRadius2_;
}
//---------------------------------------------------------------------
// return the maximum distance at which this kernel is non zero
inline double maxSquaredDistance(double h) const
{
return power<2>(getDilationFactor())*h*h;
}
private:
double norm_;
double Hinverse_;
double Hinverse2_;
double maxRadius_;
double maxRadius2_;
double scale_W_;
double scale_GradW_;
};
}}}
#endif

110
vtkm/worklet/kernels/KernelBase.h Normal file
View File

@ -0,0 +1,110 @@
#ifndef VTKM_KERNELBASE_HPP
#define VTKM_KERNELBASE_HPP
#include <vtkm/Math.h>
#include <vtkm/Types.h>
namespace vtkm { namespace worklet {
namespace kernels {
// Vector class used in the kernels
typedef vtkm::Vec<vtkm::Float64, 3> vector_type;
// templated utility to generate expansions at compile time for x^N
template <int N>
double power(double x) { return x * power<N-1>(x); }
template <>
double power<0>(double) { return 1; }
//---------------------------------------------------------------------
// Base class for Kernels
// We use CRTP to avoid virtual function calls.
template <typename Kernel>
struct KernelBase
{
//---------------------------------------------------------------------
// Constructor
// Calculate coefficients used repeatedly when evaluating the kernel
// value or gradient
// The smoothing length is usually denoted as 'h' in SPH literature
KernelBase(double smoothingLength)
: smoothingLength_(smoothingLength) {};
//---------------------------------------------------------------------
// The functions below are placeholders which should be provided by
// concrete implementations of this class.
// The KernelBase versions will not be called when algorithms are
// templated over a concrete implementation.
//---------------------------------------------------------------------
//---------------------------------------------------------------------
// compute w(h) for the given distance
inline double w(double distance) {
return static_cast<Kernel*>(this)->w(distance);
}
//---------------------------------------------------------------------
// compute w(h) for the given squared distance
// this version takes the distance squared as a convenience/optimization
// but not all implementations will benefit from it
inline double w2(double distance2) {
return static_cast<Kernel*>(this)->w2(distance2);
}
//---------------------------------------------------------------------
// compute w(h) for a variable h kernel
// this is less efficient than the fixed radius version as coefficients
// must be calculatd on the fly, but it is required when all particles
// have different smoothing lengths
inline double w(double h, double distance) {
return static_cast<Kernel*>(this)->w(h, distance);
}
//---------------------------------------------------------------------
// compute w(h) for a variable h kernel using distance squared
// this version takes the distance squared as a convenience/optimization
inline double w2(double h, double distance2) {
return static_cast<Kernel*>(this)->w2(h, distance2);
}
//---------------------------------------------------------------------
// Calculates the kernel derivative for a distance {x,y,z} vector
// from the centre.
inline vector_type gradW(double distance, const vector_type& pos) {
return static_cast<Kernel*>(this)->gradW(distance, pos);
}
// Calculates the kernel derivative at the given distance using a variable h value
// this is less efficient than the fixed radius version as coefficients
// must be calculatd on the fly
inline vector_type gradW(double h, double distance, const vector_type& pos) {
return static_cast<Kernel*>(this)->gradW(h, distance, pos);
}
// return the multiplier between smoothing length and max cutoff distance
inline double getDilationFactor() const {
return static_cast<Kernel*>(this)->getDilationFactor;
}
// return the maximum cutoff distance over which the kernel acts,
// beyond this distance the kernel value is zero
inline double maxDistance() {
return static_cast<Kernel*>(this)->maxDistance();
}
// return the maximum cutoff distance over which the kernel acts,
// beyond this distance the kernel value is zero
inline double maxDistanceSquared() {
return static_cast<Kernel*>(this)->maxDistanceSquared();
}
protected:
const double smoothingLength_;
};
}}}
#endif

199
vtkm/worklet/kernels/Spline3rdOrder.h Normal file
View File

@ -0,0 +1,199 @@
#ifndef VTKM_KERNEL_SPLINE_3RD_ORDER_H
#define VTKM_KERNEL_SPLINE_3RD_ORDER_H
#include "KernelBase.h"
//
// Spline 3rd Order kernel.
//
namespace vtkm { namespace worklet {
namespace kernels {
template <int Dimensions>
struct Spline3rdOrder : public KernelBase< Spline3rdOrder<Dimensions> >
{
//---------------------------------------------------------------------
// Constructor
// Calculate coefficients used repeatedly when evaluating the kernel
// value or gradient
Spline3rdOrder(double smoothingLength)
: KernelBase< Spline3rdOrder<Dimensions> >(smoothingLength)
{
Hinverse_ = 1.0/smoothingLength;
Hinverse2_ = Hinverse_*Hinverse_;
maxRadius_ = 2.0*smoothingLength;
maxRadius2_ = maxRadius_*maxRadius_;
//
if (Dimensions==2) {
norm_ = 10.0/(7.0*M_PI);
}
if (Dimensions==3) {
norm_ = 1.0/M_PI;
}
scale_W_ = norm_ * power<Dimensions> (Hinverse_);
scale_GradW_ = norm_ * power<Dimensions+1>(Hinverse_);
}
//---------------------------------------------------------------------
// Calculates the kernel value for the given distance
inline double w(double distance) const
{
// compute Q=(r/h)
double Q = distance * Hinverse_;
if (Q<1.0) {
return scale_W_ *(1.0 - (3.0/2.0)*Q*Q + (3.0/4.0)*Q*Q*Q);
}
else if (Q<2.0) {
double q2 = (2.0-Q);
return scale_W_ * (1.0/4.0) * (q2*q2*q2);;
}
else {
return 0.0;
}
}
//---------------------------------------------------------------------
// Calculates the kernel value for the given squared distance
inline double w2(double distance2) const
{
// compute Q
double Q = sqrt(distance2) * Hinverse_;
if (Q<1.0) {
return scale_W_ *(1.0 - (3.0/2.0)*Q*Q + (3.0/4.0)*Q*Q*Q);
}
else if (Q<2.0) {
double q2 = (2.0-Q);
return scale_W_ * (1.0/4.0) * (q2*q2*q2);;
}
else {
return 0.0;
}
}
//---------------------------------------------------------------------
// compute w(h) for a variable h kernel
inline double w(double h, double distance) const
{
double Hinverse = 1.0/h;
double scale_W = norm_ * power<Dimensions>(Hinverse);
double Q = distance * Hinverse;
if (Q<1.0) {
return scale_W *(1.0 - (3.0/2.0)*Q*Q + (3.0/4.0)*Q*Q*Q);
}
else if (Q<2.0) {
double q2 = (2.0-Q);
return scale_W * (1.0/4.0) * (q2*q2*q2);;
}
else {
return 0.0;
}
}
//---------------------------------------------------------------------
// compute w(h) for a variable h kernel using distance squared
inline double w2(double h, double distance2) const
{
double Hinverse = 1.0/h;
double scale_W = norm_ * power<Dimensions>(Hinverse);
double Q = sqrt(distance2) * Hinverse;
if (Q<1.0) {
return scale_W *(1.0 - (3.0/2.0)*Q*Q + (3.0/4.0)*Q*Q*Q);
}
else if (Q<2.0) {
double q2 = (2.0-Q);
return scale_W * (1.0/4.0) * (q2*q2*q2);;
}
else {
return 0.0;
}
}
//---------------------------------------------------------------------
// Calculates the kernel derivation for the given distance of two particles.
// The used formula is the derivation of Speith (3.126) for the value
// with (3.21) for the direction of the gradient vector.
// Be careful: grad W is antisymmetric in r (3.25)!.
inline vector_type gradW(double distance, const vector_type& pos) const
{
double Q = distance * Hinverse_;
if (Q==0.0) {
return vector_type(0.0);
}
else if (Q<1.0) {
return scale_GradW_ * (-3.0*Q + (9.0/4.0)*Q*Q) * pos;
}
else if (Q<2.0) {
double q2 = (2.0-Q);
return scale_GradW_ * (-3.0/4.0)*q2*q2 * pos;
}
else {
return vector_type(0.0);
}
}
//---------------------------------------------------------------------
inline vector_type gradW(double h, double distance, const vector_type& pos) const
{
double Hinverse = 1.0/h;
double scale_GradW = norm_ * power<Dimensions+1>(Hinverse);
double Q = distance * Hinverse;
if (Q==0.0) {
return vector_type(0.0);
}
else if (Q<1.0) {
return scale_GradW * (-3.0*Q + (9.0/4.0)*Q*Q) * pos;
}
else if (Q<2.0) {
double q2 = (2.0-Q);
return scale_GradW * (-3.0/4.0)*q2*q2 * pos;
}
else {
return vector_type(0.0);
}
}
//---------------------------------------------------------------------
// return the maximum distance at which this kernel is non zero
inline double maxDistance() const
{
return maxRadius_;
}
//---------------------------------------------------------------------
// return the maximum distance at which this variable h kernel is non zero
inline double maxDistance(double h) const
{
return 2.0*h;
}
//---------------------------------------------------------------------
// return the maximum distance at which this kernel is non zero
inline double maxSquaredDistance() const
{
return maxRadius2_;
}
//---------------------------------------------------------------------
// return the maximum distance at which this kernel is non zero
inline double maxSquaredDistance(double h) const
{
return 4.0*h*h;
}
//---------------------------------------------------------------------
// return the multiplier between smoothing length and max cutoff distance
inline double getDilationFactor() const { return 2.0; }
private:
double norm_;
double Hinverse_;
double Hinverse2_;
double maxRadius_;
double maxRadius2_;
double scale_W_;
double scale_GradW_;
};
}}}
#endif