You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
260 lines
8.4 KiB
260 lines
8.4 KiB
/*=========================================================================
|
|
|
|
Program: Visualization Toolkit
|
|
Module: $RCSfile: vtkPointLoad.cxx,v $
|
|
|
|
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
|
|
All rights reserved.
|
|
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
|
|
|
|
This software is distributed WITHOUT ANY WARRANTY; without even
|
|
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
|
|
PURPOSE. See the above copyright notice for more information.
|
|
|
|
=========================================================================*/
|
|
#include "vtkPointLoad.h"
|
|
|
|
#include "vtkFloatArray.h"
|
|
#include "vtkImageData.h"
|
|
#include "vtkMath.h"
|
|
#include "vtkInformation.h"
|
|
#include "vtkInformationVector.h"
|
|
#include "vtkObjectFactory.h"
|
|
#include "vtkStreamingDemandDrivenPipeline.h"
|
|
#include "vtkPointData.h"
|
|
|
|
vtkCxxRevisionMacro(vtkPointLoad, "$Revision: 1.49 $");
|
|
vtkStandardNewMacro(vtkPointLoad);
|
|
|
|
// Construct with ModelBounds=(-1,1,-1,1,-1,1), SampleDimensions=(50,50,50),
|
|
// and LoadValue = 1.
|
|
vtkPointLoad::vtkPointLoad()
|
|
{
|
|
this->LoadValue = 1.0;
|
|
|
|
this->ModelBounds[0] = -1.0;
|
|
this->ModelBounds[1] = 1.0;
|
|
this->ModelBounds[2] = -1.0;
|
|
this->ModelBounds[3] = 1.0;
|
|
this->ModelBounds[4] = -1.0;
|
|
this->ModelBounds[5] = 1.0;
|
|
|
|
this->SampleDimensions[0] = 50;
|
|
this->SampleDimensions[1] = 50;
|
|
this->SampleDimensions[2] = 50;
|
|
|
|
this->PoissonsRatio = 0.3;
|
|
|
|
this->SetNumberOfInputPorts(0);
|
|
}
|
|
|
|
// Specify the dimensions of the volume. A stress tensor will be computed for
|
|
// each point in the volume.
|
|
void vtkPointLoad::SetSampleDimensions(int i, int j, int k)
|
|
{
|
|
int dim[3];
|
|
|
|
dim[0] = i;
|
|
dim[1] = j;
|
|
dim[2] = k;
|
|
|
|
this->SetSampleDimensions(dim);
|
|
}
|
|
|
|
// Specify the dimensions of the volume. A stress tensor will be computed for
|
|
// each point in the volume.
|
|
void vtkPointLoad::SetSampleDimensions(int dim[3])
|
|
{
|
|
vtkDebugMacro(<< " setting SampleDimensions to (" << dim[0] << "," << dim[1] << "," << dim[2] << ")");
|
|
|
|
if ( dim[0] != this->SampleDimensions[0] ||
|
|
dim[1] != this->SampleDimensions[1] ||
|
|
dim[2] != this->SampleDimensions[2] )
|
|
{
|
|
for ( int i=0; i<3; i++)
|
|
{
|
|
this->SampleDimensions[i] = (dim[i] > 0 ? dim[i] : 1);
|
|
}
|
|
this->Modified();
|
|
}
|
|
}
|
|
|
|
int vtkPointLoad::RequestInformation (
|
|
vtkInformation * vtkNotUsed(request),
|
|
vtkInformationVector ** vtkNotUsed( inputVector ),
|
|
vtkInformationVector *outputVector)
|
|
{
|
|
// get the info objects
|
|
vtkInformation* outInfo = outputVector->GetInformationObject(0);
|
|
|
|
// use model bounds
|
|
double origin[3];
|
|
origin[0] = this->ModelBounds[0];
|
|
origin[1] = this->ModelBounds[2];
|
|
origin[2] = this->ModelBounds[4];
|
|
outInfo->Set(vtkDataObject::ORIGIN(), origin, 3);
|
|
|
|
// Set volume origin and data spacing
|
|
int i;
|
|
double spacing[3];
|
|
for (i=0; i<3; i++)
|
|
{
|
|
spacing[i] = (this->ModelBounds[2*i+1] - this->ModelBounds[2*i])
|
|
/ (this->SampleDimensions[i] - 1);
|
|
if ( spacing[i] <= 0.0 )
|
|
{
|
|
spacing[i] = 1.0;
|
|
}
|
|
}
|
|
outInfo->Set(vtkDataObject::SPACING(),spacing,3);
|
|
|
|
int wExt[6];
|
|
wExt[0] = 0; wExt[2] = 0; wExt[4] = 0;
|
|
wExt[1] = this->SampleDimensions[0] - 1;
|
|
wExt[3] = this->SampleDimensions[1] - 1;
|
|
wExt[5] = this->SampleDimensions[2] - 1;
|
|
|
|
outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), wExt, 6);
|
|
vtkDataObject::SetPointDataActiveScalarInfo(outInfo, VTK_FLOAT, 1);
|
|
return 1;
|
|
}
|
|
|
|
//
|
|
// Generate tensors and scalars for point load on semi-infinite domain.
|
|
//
|
|
void vtkPointLoad::ExecuteData(vtkDataObject *outp)
|
|
{
|
|
int i, j, k;
|
|
vtkFloatArray *newTensors;
|
|
double tensor[9];
|
|
vtkIdType numPts;
|
|
double P, twoPi, xP[3], rho, rho2, rho3, rho5, nu;
|
|
double x, x2, y, y2, z, z2, rhoPlusz2, zPlus2rho, txy, txz, tyz;
|
|
double sx, sy, sz, seff;
|
|
vtkImageData *output = this->AllocateOutputData(outp);
|
|
vtkFloatArray *newScalars =
|
|
vtkFloatArray::SafeDownCast(output->GetPointData()->GetScalars());
|
|
double *spacing, *origin;
|
|
|
|
vtkDebugMacro(<< "Computing point load stress tensors");
|
|
|
|
//
|
|
// Initialize self; create output objects
|
|
//
|
|
numPts = this->SampleDimensions[0] * this->SampleDimensions[1]
|
|
* this->SampleDimensions[2];
|
|
spacing = output->GetSpacing();
|
|
origin = output->GetOrigin();
|
|
newTensors = vtkFloatArray::New();
|
|
newTensors->SetNumberOfComponents(9);
|
|
newTensors->Allocate(9*numPts);
|
|
|
|
//
|
|
// Compute the location of the load
|
|
//
|
|
xP[0] = (this->ModelBounds[0] + this->ModelBounds[1]) / 2.0; //in center
|
|
xP[1] = (this->ModelBounds[2] + this->ModelBounds[3]) / 2.0;
|
|
xP[2] = this->ModelBounds[5]; // at top of box
|
|
//
|
|
// Traverse all points evaluating implicit function at each point. Note that
|
|
// points are evaluated in local coordinate system of applied force.
|
|
//
|
|
twoPi = 2.0*vtkMath::Pi();
|
|
P = -this->LoadValue;
|
|
int pointCount = 0;
|
|
for (k=0; k<this->SampleDimensions[2]; k++)
|
|
{
|
|
z = xP[2] - (origin[2] + k*spacing[2]);
|
|
for (j=0; j<this->SampleDimensions[1]; j++)
|
|
{
|
|
y = xP[1] - (origin[1] + j*spacing[1]);
|
|
for (i=0; i<this->SampleDimensions[0]; i++)
|
|
{
|
|
x = (origin[0] + i*spacing[0]) - xP[0];
|
|
rho = sqrt(x*x + y*y + z*z);//in local coordinates
|
|
if ( rho < 1.0e-10 )
|
|
{
|
|
vtkWarningMacro(<<"Attempting to set singularity, resetting");
|
|
tensor[0] = VTK_LARGE_FLOAT; // Component(0,0)
|
|
tensor[4] = VTK_LARGE_FLOAT; // Component(1,1);
|
|
tensor[8] = VTK_LARGE_FLOAT; // Component(2,2);
|
|
tensor[3] = 0.0; // Component(0,1);
|
|
tensor[6] = 0.0; // Component(0,2);
|
|
tensor[1] = 0.0; // Component(1,0);
|
|
tensor[7] = 0.0; // Component(1,2);
|
|
tensor[2] = 0.0; // Component(2,0);
|
|
tensor[5] = 0.0; // Component(2,1);
|
|
newTensors->InsertNextTuple(tensor);
|
|
double val = VTK_LARGE_FLOAT;
|
|
newScalars->InsertTuple(pointCount,&val);
|
|
pointCount++;
|
|
continue;
|
|
}
|
|
|
|
rho2 = rho*rho;
|
|
rho3 = rho2*rho;
|
|
rho5 = rho2*rho3;
|
|
nu = (1.0 - 2.0*this->PoissonsRatio);
|
|
x2 = x*x;
|
|
y2 = y*y;
|
|
z2 = z*z;
|
|
rhoPlusz2 = (rho + z) * (rho + z);
|
|
zPlus2rho = (2.0*rho + z);
|
|
|
|
// normal stresses
|
|
sx = P/(twoPi*rho2) * (3.0*z*x2/rho3 - nu*(z/rho - rho/(rho+z) +
|
|
x2*(zPlus2rho)/(rho*rhoPlusz2)));
|
|
sy = P/(twoPi*rho2) * (3.0*z*y2/rho3 - nu*(z/rho - rho/(rho+z) +
|
|
y2*(zPlus2rho)/(rho*rhoPlusz2)));
|
|
sz = 3.0*P*z2*z/(twoPi*rho5);
|
|
|
|
//shear stresses - negative signs are coordinate transformations
|
|
//that is, equations (in text) are in different coordinate system
|
|
//than volume is in.
|
|
txy = -(P/(twoPi*rho2) * (3.0*x*y*z/rho3 -
|
|
nu*x*y*(zPlus2rho)/(rho*rhoPlusz2)));
|
|
txz = -(3.0*P*x*z2/(twoPi*rho5));
|
|
tyz = 3.0*P*y*z2/(twoPi*rho5);
|
|
|
|
tensor[0] = sx; // Component(0,0);
|
|
tensor[4] = sy; // Component(1,1);
|
|
tensor[8] = sz; // Component(2,2);
|
|
tensor[3] = txy; // Component(0,1); real symmetric matrix
|
|
tensor[1] = txy; // Component(1,0);
|
|
tensor[6] = txz; // Component(0,2);
|
|
tensor[2] = txz; // Component(2,0);
|
|
tensor[7] = tyz; // Component(1,2);
|
|
tensor[5] = tyz; // Component(2,1);
|
|
newTensors->InsertNextTuple(tensor);
|
|
|
|
seff = 0.333333* sqrt ((sx-sy)*(sx-sy) + (sy-sz)*(sy-sz) +
|
|
(sz-sx)*(sz-sx) + 6.0*txy*txy + 6.0*tyz*tyz +
|
|
6.0*txz*txz);
|
|
newScalars->InsertTuple(pointCount,&seff);
|
|
pointCount++;
|
|
}
|
|
}
|
|
}
|
|
//
|
|
// Update self and free memory
|
|
//
|
|
output->GetPointData()->SetTensors(newTensors);
|
|
newTensors->Delete();
|
|
}
|
|
|
|
|
|
void vtkPointLoad::PrintSelf(ostream& os, vtkIndent indent)
|
|
{
|
|
this->Superclass::PrintSelf(os,indent);
|
|
|
|
os << indent << "Load Value: " << this->LoadValue << "\n";
|
|
os << indent << "Sample Dimensions: (" << this->SampleDimensions[0] << ", "
|
|
<< this->SampleDimensions[1] << ", "
|
|
<< this->SampleDimensions[2] << ")\n";
|
|
os << indent << "ModelBounds: \n";
|
|
os << indent << " Xmin,Xmax: (" << this->ModelBounds[0] << ", " << this->ModelBounds[1] << ")\n";
|
|
os << indent << " Ymin,Ymax: (" << this->ModelBounds[2] << ", " << this->ModelBounds[3] << ")\n";
|
|
os << indent << " Zmin,Zmax: (" << this->ModelBounds[4] << ", " << this->ModelBounds[5] << ")\n";
|
|
os << indent << "Poisson's Ratio: " << this->PoissonsRatio << "\n";
|
|
}
|
|
|
|
|