1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkPointLoad.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 "vtkPointLoad.h"
16
17 #include "vtkFloatArray.h"
18 #include "vtkImageData.h"
19 #include "vtkMath.h"
20 #include "vtkInformation.h"
21 #include "vtkInformationVector.h"
22 #include "vtkObjectFactory.h"
23 #include "vtkStreamingDemandDrivenPipeline.h"
24 #include "vtkPointData.h"
25
26 vtkStandardNewMacro(vtkPointLoad);
27
28 // Construct with ModelBounds=(-1,1,-1,1,-1,1), SampleDimensions=(50,50,50),
29 // and LoadValue = 1.
vtkPointLoad()30 vtkPointLoad::vtkPointLoad()
31 {
32 this->LoadValue = 1.0;
33
34 this->ModelBounds[0] = -1.0;
35 this->ModelBounds[1] = 1.0;
36 this->ModelBounds[2] = -1.0;
37 this->ModelBounds[3] = 1.0;
38 this->ModelBounds[4] = -1.0;
39 this->ModelBounds[5] = 1.0;
40
41 this->SampleDimensions[0] = 50;
42 this->SampleDimensions[1] = 50;
43 this->SampleDimensions[2] = 50;
44
45 this->PoissonsRatio = 0.3;
46
47 this->SetNumberOfInputPorts(0);
48 }
49
50 // Specify the dimensions of the volume. A stress tensor will be computed for
51 // each point in the volume.
SetSampleDimensions(int i,int j,int k)52 void vtkPointLoad::SetSampleDimensions(int i, int j, int k)
53 {
54 int dim[3];
55
56 dim[0] = i;
57 dim[1] = j;
58 dim[2] = k;
59
60 this->SetSampleDimensions(dim);
61 }
62
63 // Specify the dimensions of the volume. A stress tensor will be computed for
64 // each point in the volume.
SetSampleDimensions(int dim[3])65 void vtkPointLoad::SetSampleDimensions(int dim[3])
66 {
67 vtkDebugMacro(<< " setting SampleDimensions to (" << dim[0] << "," << dim[1] << "," << dim[2] << ")");
68
69 if ( dim[0] != this->SampleDimensions[0] ||
70 dim[1] != this->SampleDimensions[1] ||
71 dim[2] != this->SampleDimensions[2] )
72 {
73 for ( int i=0; i<3; i++)
74 {
75 this->SampleDimensions[i] = (dim[i] > 0 ? dim[i] : 1);
76 }
77 this->Modified();
78 }
79 }
80
RequestInformation(vtkInformation * vtkNotUsed (request),vtkInformationVector ** vtkNotUsed (inputVector),vtkInformationVector * outputVector)81 int vtkPointLoad::RequestInformation (
82 vtkInformation * vtkNotUsed(request),
83 vtkInformationVector ** vtkNotUsed( inputVector ),
84 vtkInformationVector *outputVector)
85 {
86 // get the info objects
87 vtkInformation* outInfo = outputVector->GetInformationObject(0);
88
89 // use model bounds
90 double origin[3];
91 origin[0] = this->ModelBounds[0];
92 origin[1] = this->ModelBounds[2];
93 origin[2] = this->ModelBounds[4];
94 outInfo->Set(vtkDataObject::ORIGIN(), origin, 3);
95
96 // Set volume origin and data spacing
97 int i;
98 double spacing[3];
99 for (i=0; i<3; i++)
100 {
101 spacing[i] = (this->ModelBounds[2*i+1] - this->ModelBounds[2*i])
102 / (this->SampleDimensions[i] - 1);
103 if ( spacing[i] <= 0.0 )
104 {
105 spacing[i] = 1.0;
106 }
107 }
108 outInfo->Set(vtkDataObject::SPACING(),spacing,3);
109
110 int wExt[6];
111 wExt[0] = 0; wExt[2] = 0; wExt[4] = 0;
112 wExt[1] = this->SampleDimensions[0] - 1;
113 wExt[3] = this->SampleDimensions[1] - 1;
114 wExt[5] = this->SampleDimensions[2] - 1;
115
116 outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), wExt, 6);
117 vtkDataObject::SetPointDataActiveScalarInfo(outInfo, VTK_FLOAT, 1);
118 return 1;
119 }
120
121 //
122 // Generate tensors and scalars for point load on semi-infinite domain.
123 //
ExecuteDataWithInformation(vtkDataObject * outp,vtkInformation * outInfo)124 void vtkPointLoad::ExecuteDataWithInformation(vtkDataObject *outp, vtkInformation* outInfo)
125 {
126 int i, j, k;
127 vtkFloatArray *newTensors;
128 double tensor[9];
129 vtkIdType numPts;
130 double P, twoPi, xP[3], rho, rho2, rho3, rho5, nu;
131 double x, x2, y, y2, z, z2, rhoPlusz2, zPlus2rho, txy, txz, tyz;
132 double sx, sy, sz, seff;
133 vtkImageData *output = this->AllocateOutputData(outp, outInfo);
134 vtkFloatArray *newScalars =
135 vtkFloatArray::SafeDownCast(output->GetPointData()->GetScalars());
136 double *spacing, *origin;
137
138 vtkDebugMacro(<< "Computing point load stress tensors");
139
140 //
141 // Initialize self; create output objects
142 //
143 numPts = this->SampleDimensions[0] * this->SampleDimensions[1]
144 * this->SampleDimensions[2];
145 spacing = output->GetSpacing();
146 origin = output->GetOrigin();
147 newTensors = vtkFloatArray::New();
148 newTensors->SetNumberOfComponents(9);
149 newTensors->Allocate(9*numPts);
150
151 //
152 // Compute the location of the load
153 //
154 xP[0] = (this->ModelBounds[0] + this->ModelBounds[1]) / 2.0; //in center
155 xP[1] = (this->ModelBounds[2] + this->ModelBounds[3]) / 2.0;
156 xP[2] = this->ModelBounds[5]; // at top of box
157 //
158 // Traverse all points evaluating implicit function at each point. Note that
159 // points are evaluated in local coordinate system of applied force.
160 //
161 twoPi = 2.0*vtkMath::Pi();
162 P = -this->LoadValue;
163 int pointCount = 0;
164 for (k=0; k<this->SampleDimensions[2]; k++)
165 {
166 z = xP[2] - (origin[2] + k*spacing[2]);
167 for (j=0; j<this->SampleDimensions[1]; j++)
168 {
169 y = xP[1] - (origin[1] + j*spacing[1]);
170 for (i=0; i<this->SampleDimensions[0]; i++)
171 {
172 x = (origin[0] + i*spacing[0]) - xP[0];
173 rho = sqrt(x*x + y*y + z*z);//in local coordinates
174 if ( rho < 1.0e-10 )
175 {
176 vtkWarningMacro(<<"Attempting to set singularity, resetting");
177 tensor[0] = VTK_FLOAT_MAX; // Component(0,0)
178 tensor[4] = VTK_FLOAT_MAX; // Component(1,1);
179 tensor[8] = VTK_FLOAT_MAX; // Component(2,2);
180 tensor[3] = 0.0; // Component(0,1);
181 tensor[6] = 0.0; // Component(0,2);
182 tensor[1] = 0.0; // Component(1,0);
183 tensor[7] = 0.0; // Component(1,2);
184 tensor[2] = 0.0; // Component(2,0);
185 tensor[5] = 0.0; // Component(2,1);
186 newTensors->InsertNextTuple(tensor);
187 double val = VTK_FLOAT_MAX;
188 newScalars->InsertTuple(pointCount,&val);
189 pointCount++;
190 continue;
191 }
192
193 rho2 = rho*rho;
194 rho3 = rho2*rho;
195 rho5 = rho2*rho3;
196 nu = (1.0 - 2.0*this->PoissonsRatio);
197 x2 = x*x;
198 y2 = y*y;
199 z2 = z*z;
200 rhoPlusz2 = (rho + z) * (rho + z);
201 zPlus2rho = (2.0*rho + z);
202
203 // normal stresses
204 sx = P/(twoPi*rho2) * (3.0*z*x2/rho3 - nu*(z/rho - rho/(rho+z) +
205 x2*(zPlus2rho)/(rho*rhoPlusz2)));
206 sy = P/(twoPi*rho2) * (3.0*z*y2/rho3 - nu*(z/rho - rho/(rho+z) +
207 y2*(zPlus2rho)/(rho*rhoPlusz2)));
208 sz = 3.0*P*z2*z/(twoPi*rho5);
209
210 //shear stresses - negative signs are coordinate transformations
211 //that is, equations (in text) are in different coordinate system
212 //than volume is in.
213 txy = -(P/(twoPi*rho2) * (3.0*x*y*z/rho3 -
214 nu*x*y*(zPlus2rho)/(rho*rhoPlusz2)));
215 txz = -(3.0*P*x*z2/(twoPi*rho5));
216 tyz = 3.0*P*y*z2/(twoPi*rho5);
217
218 tensor[0] = sx; // Component(0,0);
219 tensor[4] = sy; // Component(1,1);
220 tensor[8] = sz; // Component(2,2);
221 tensor[3] = txy; // Component(0,1); real symmetric matrix
222 tensor[1] = txy; // Component(1,0);
223 tensor[6] = txz; // Component(0,2);
224 tensor[2] = txz; // Component(2,0);
225 tensor[7] = tyz; // Component(1,2);
226 tensor[5] = tyz; // Component(2,1);
227 newTensors->InsertNextTuple(tensor);
228
229 seff = 0.333333* sqrt ((sx-sy)*(sx-sy) + (sy-sz)*(sy-sz) +
230 (sz-sx)*(sz-sx) + 6.0*txy*txy + 6.0*tyz*tyz +
231 6.0*txz*txz);
232 newScalars->InsertTuple(pointCount,&seff);
233 pointCount++;
234 }
235 }
236 }
237 //
238 // Update self and free memory
239 //
240 output->GetPointData()->SetTensors(newTensors);
241 newTensors->Delete();
242 }
243
244
PrintSelf(ostream & os,vtkIndent indent)245 void vtkPointLoad::PrintSelf(ostream& os, vtkIndent indent)
246 {
247 this->Superclass::PrintSelf(os,indent);
248
249 os << indent << "Load Value: " << this->LoadValue << "\n";
250 os << indent << "Sample Dimensions: (" << this->SampleDimensions[0] << ", "
251 << this->SampleDimensions[1] << ", "
252 << this->SampleDimensions[2] << ")\n";
253 os << indent << "ModelBounds: \n";
254 os << indent << " Xmin,Xmax: (" << this->ModelBounds[0] << ", " << this->ModelBounds[1] << ")\n";
255 os << indent << " Ymin,Ymax: (" << this->ModelBounds[2] << ", " << this->ModelBounds[3] << ")\n";
256 os << indent << " Zmin,Zmax: (" << this->ModelBounds[4] << ", " << this->ModelBounds[5] << ")\n";
257 os << indent << "Poisson's Ratio: " << this->PoissonsRatio << "\n";
258 }
259
260