1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkSPHKernel.cxx
5
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 All rights reserved.
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
13
14 =========================================================================*/
15 #include "vtkSPHKernel.h"
16 #include "vtkAbstractPointLocator.h"
17 #include "vtkObjectFactory.h"
18 #include "vtkIdList.h"
19 #include "vtkDoubleArray.h"
20 #include "vtkFloatArray.h"
21 #include "vtkDataArray.h"
22 #include "vtkDataSet.h"
23 #include "vtkMath.h"
24
25 vtkCxxSetObjectMacro(vtkSPHKernel,CutoffArray,vtkDataArray);
26 vtkCxxSetObjectMacro(vtkSPHKernel,DensityArray,vtkDataArray);
27 vtkCxxSetObjectMacro(vtkSPHKernel,MassArray,vtkDataArray);
28
29 //----------------------------------------------------------------------------
vtkSPHKernel()30 vtkSPHKernel::vtkSPHKernel()
31 {
32 this->RequiresInitialization = true;
33 this->SpatialStep = 0.001;
34 this->Dimension = 3;
35 this->CutoffArray = nullptr;
36 this->DensityArray = nullptr;
37 this->MassArray = nullptr;
38 }
39
40 //----------------------------------------------------------------------------
~vtkSPHKernel()41 vtkSPHKernel::~vtkSPHKernel()
42 {
43 this->SetCutoffArray(nullptr);
44 this->SetDensityArray(nullptr);
45 this->SetMassArray(nullptr);
46 }
47
48 //----------------------------------------------------------------------------
49 // At this point, the spatial step, the dimension of the kernel, the cutoff
50 // factor, and the sigma normalization factor should be known.
51 void vtkSPHKernel::
Initialize(vtkAbstractPointLocator * loc,vtkDataSet * ds,vtkPointData * attr)52 Initialize(vtkAbstractPointLocator *loc, vtkDataSet *ds, vtkPointData *attr)
53 {
54 this->Superclass::Initialize(loc, ds, attr);
55
56 // this->CutoffFactor should have been set by subclass
57 this->Cutoff = this->CutoffFactor * this->SpatialStep;
58 this->DistNorm = 1.0 / this->SpatialStep;
59 this->NormFactor = this->Sigma * pow(this->DistNorm,this->Dimension);
60 this->DefaultVolume = pow(this->SpatialStep,this->Dimension);
61
62 // See if cutoff array is provided.
63 if ( this->CutoffArray && this->CutoffArray->GetNumberOfComponents() == 1 )
64 {
65 this->UseCutoffArray = true;
66 }
67 else
68 {
69 this->UseCutoffArray = false;
70 }
71
72 // See if local mass and density information is provided
73 if ( this->DensityArray && this->MassArray &&
74 this->DensityArray->GetNumberOfComponents() == 1 &&
75 this->MassArray->GetNumberOfComponents() == 1 )
76 {
77 this->UseArraysForVolume = true;
78 }
79 else
80 {
81 this->UseArraysForVolume = false;
82 }
83 }
84
85 //----------------------------------------------------------------------------
86 // Radius around point is cutoff factor * smoothing length. That is unless
87 // cutoff array is provided.
88 vtkIdType vtkSPHKernel::
ComputeBasis(double x[3],vtkIdList * pIds,vtkIdType ptId)89 ComputeBasis(double x[3], vtkIdList *pIds, vtkIdType ptId)
90 {
91 double cutoff;
92 if ( this->UseCutoffArray )
93 {
94 this->CutoffArray->GetTuple(ptId,&cutoff);
95 }
96 else
97 {
98 cutoff = this->Cutoff;
99 }
100
101 this->Locator->FindPointsWithinRadius(cutoff, x, pIds);
102 return pIds->GetNumberOfIds();
103 }
104
105 //----------------------------------------------------------------------------
106 vtkIdType vtkSPHKernel::
ComputeWeights(double x[3],vtkIdList * pIds,vtkDoubleArray * weights)107 ComputeWeights(double x[3], vtkIdList *pIds, vtkDoubleArray *weights)
108 {
109 vtkIdType numPts = pIds->GetNumberOfIds();
110 int i;
111 vtkIdType id;
112 double d, y[3];
113 weights->SetNumberOfTuples(numPts);
114 double *w = weights->GetPointer(0);
115 double KW, mass, density, volume;
116
117 // Compute SPH coefficients.
118 for (i=0; i<numPts; ++i)
119 {
120 id = pIds->GetId(i);
121 this->DataSet->GetPoint(id,y);
122 d = sqrt( vtkMath::Distance2BetweenPoints(x,y) );
123
124 KW = this->ComputeFunctionWeight(d*this->DistNorm);
125
126 if ( this->UseArraysForVolume )
127 {
128 this->MassArray->GetTuple(id,&mass);
129 this->DensityArray->GetTuple(id,&density);
130 volume = mass /density;
131 }
132 else
133 {
134 volume = this->DefaultVolume;
135 }
136
137 w[i] = this->NormFactor * KW * volume;
138 }//over all neighbor points
139
140 return numPts;
141 }
142
143 //----------------------------------------------------------------------------
144 vtkIdType vtkSPHKernel::
ComputeDerivWeights(double x[3],vtkIdList * pIds,vtkDoubleArray * weights,vtkDoubleArray * gradWeights)145 ComputeDerivWeights(double x[3], vtkIdList *pIds, vtkDoubleArray *weights,
146 vtkDoubleArray *gradWeights)
147 {
148 vtkIdType numPts = pIds->GetNumberOfIds();
149 int i;
150 vtkIdType id;
151 double d, y[3];
152 weights->SetNumberOfTuples(numPts);
153 double *w = weights->GetPointer(0);
154 gradWeights->SetNumberOfTuples(numPts);
155 double *gw = gradWeights->GetPointer(0);
156 double KW, GW, volume=this->DefaultVolume;
157
158 // Compute SPH coefficients for data and deriative data
159 for (i=0; i<numPts; ++i)
160 {
161 id = pIds->GetId(i);
162 this->DataSet->GetPoint(id,y);
163 d = sqrt( vtkMath::Distance2BetweenPoints(x,y) );
164
165 KW = this->ComputeFunctionWeight(d*this->DistNorm);
166 GW = this->ComputeDerivWeight(d*this->DistNorm);
167
168 w[i] = this->NormFactor * KW * volume;
169 gw[i] = this->NormFactor * GW * volume;
170 }//over all neighbor points
171
172 return numPts;
173 }
174
175 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)176 void vtkSPHKernel::PrintSelf(ostream& os, vtkIndent indent)
177 {
178 this->Superclass::PrintSelf(os,indent);
179
180 os << indent << "Spatial Step: " << this->SpatialStep << "\n";
181 os << indent << "Dimension: " << this->Dimension << "\n";
182 os << indent << "Cutoff Factor: " << this->CutoffFactor << "\n";
183 os << indent << "Sigma: " << this->Sigma << "\n";
184
185 os << indent << "Cutoff Array: " << this->CutoffArray<< "\n";
186 os << indent << "Density Array: " << this->DensityArray << "\n";
187 os << indent << "Mass Array: " << this->MassArray << "\n";
188 }
189