Commit 62d9cce4 authored by Will Schroeder's avatar Will Schroeder
Browse files

Initial revision

parent 3c4c1925
/*=========================================================================
Program: Visualization Library
Module: CellType.hh
Language: C++
Date: $Date$
Version: $Revision$
Description:
---------------------------------------------------------------------------
This file is part of the Visualization Library. No part of this file
or its contents may be copied, reproduced or altered in any way
without the express written consent of the authors.
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen 1993, 1994
=========================================================================*/
//
// Class to localize mapping of cell types to object. Adding new cell type
// means adding new & unique #define and then implementing the object.
//
#ifndef __vlCellTypes_h
#define __vlCellTypes_h
class vlCell;
#define vlPOINT 0
#define vlPOLY_POINTS 1
#define vlLINE 2
#define vlPOLY_LINE 3
#define vlTRIANGLE 4
#define vlTRIANGLE_STRIP 5
#define vlPOLYGON 6
#define vlRECTANGLE 7
#define vlQUAD 8
#define vlTETRA 9
#define vlBRICK 10
#define vlHEXAHEDRON 11
class vlCellTypes
{
public:
vlCell *MapType(int type);
};
#endif
/*=========================================================================
Program: Visualization Library
Module: Hexa.hh
Language: C++
Date: $Date$
Version: $Revision$
Description:
---------------------------------------------------------------------------
This file is part of the Visualization Library. No part of this file
or its contents may be copied, reproduced or altered in any way
without the express written consent of the authors.
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen 1993, 1994
=========================================================================*/
//
// Computational class for hexahedron
//
#ifndef __vlHexahedron_h
#define __vlHexahedron_h
#include "Cell.hh"
class vlHexahedron : public vlCell
{
public:
vlHexahedron() {};
char *GetClassName() {return "vlHexahedron";};
int CellDimension() {return 3;};
float EvaluatePosition(float x[3], int& subId, float pcoords[3]);
void EvaluateLocation(int& subId, float pcoords[3], float x[3]);
void ShapeFunctions(float pcoords[3], float sf[8]);
void ShapeDerivs(float pcoords[3], float derivs[24]);
};
#endif
/*=========================================================================
Program: Visualization Library
Module: Pixel.hh
Language: C++
Date: $Date$
Version: $Revision$
Description:
---------------------------------------------------------------------------
This file is part of the Visualization Library. No part of this file
or its contents may be copied, reproduced or altered in any way
without the express written consent of the authors.
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen 1993, 1994
=========================================================================*/
//
// Computational class for rectangles.
//
#ifndef __vlRectangle_h
#define __vlRectangle_h
#include "Cell.hh"
class vlRectangle : public vlCell
{
public:
vlRectangle() {};
char *GetClassName() {return "vlRectangle";};
int CellDimension() {return 2;};
float EvaluatePosition(float x[3], int& subId, float pcoords[3]);
void EvaluateLocation(int& subId, float pcoords[3], float x[3]);
};
#endif
/*=========================================================================
Program: Visualization Library
Module: CellArr.cc
Language: C++
Date: $Date$
Version: $Revision$
Description:
---------------------------------------------------------------------------
This file is part of the Visualization Library. No part of this file
or its contents may be copied, reproduced or altered in any way
without the express written consent of the authors.
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen 1993, 1994
=========================================================================*/
#include "CellArr.hh"
vlCellArray::vlCellArray (const vlCellArray& ca)
{
this->NumberOfCells = ca.NumberOfCells;
this->Location = 0;
this->Ia = ca.Ia;
}
/*=========================================================================
Program: Visualization Library
Module: Hexa.cc
Language: C++
Date: $Date$
Version: $Revision$
Description:
---------------------------------------------------------------------------
This file is part of the Visualization Library. No part of this file
or its contents may be copied, reproduced or altered in any way
without the express written consent of the authors.
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen 1993, 1994
=========================================================================*/
#include <math.h>
#include "Hexa.hh"
#include "vlMath.hh"
//
// Note: the ordering of the Points and PointIds is important. See text.
//
//
// Method to calculate parametric coordinates in an eight noded
// linear hexahedron element from global coordinates. Note: the natural
// formulation calls for r,s,t parametric coordinates to range range
// from -1<=r,s,t<=1. We need to shift to 0<=r,s,t<=1.
// Uses Newton's method to solve the non-linear equations.
//
#define MAX_ITERATION 10
#define CONVERGED 1.e-03
float vlHexahedron::EvaluatePosition(float x[3], int& subId, float pcoords[3])
{
int iteration, converged;
float params[3];
float fcol[3], rcol[3], scol[3], tcol[3];
int i,j,numPts,idx;
float d, *pt;
float closestPoint[3];
vlMath math;
float sf[8], derivs[24];
//
// set initial position for Newton's method
//
subId = 0;
pcoords[0] = pcoords[1] = pcoords[2] = params[0] = params[1] = params[2] = 0.0;
//
// enter iteration loop
///
for (iteration=converged=0; !converged && (iteration < MAX_ITERATION);
iteration++)
{
//
// calculate element shape functions and derivatives
//
this->ShapeFunctions(pcoords, sf);
this->ShapeDerivs(pcoords, derivs);
//
// calculate newton functions
//
for (i=0; i<3; i++)
{
fcol[i] = rcol[i] = scol[i] = tcol[i] = 0.0;
}
for (i=0; i<8; i++)
{
pt = this->Points.GetPoint(i);
for (j=0; j<3; j++)
{
fcol[j] += pt[j] * sf[i];
rcol[j] += pt[j] * derivs[i];
scol[j] += pt[j] * derivs[i+8];
tcol[j] += pt[j] * derivs[i+16];
}
}
for (i=0; i<3; i++) fcol[i] -= x[i];
//
// compute determinates and generate improvements
//
if ( (d=math.Determinate3x3(rcol,scol,tcol)) == 0.0 )
{
return LARGE_FLOAT;
}
pcoords[0] = params[0] - math.Determinate3x3 (fcol,scol,tcol) / d;
pcoords[1] = params[1] - math.Determinate3x3 (rcol,fcol,tcol) / d;
pcoords[2] = params[2] - math.Determinate3x3 (rcol,scol,fcol) / d;
//
// check for convergence
//
if ( ((fabs(pcoords[0]-params[0])) < CONVERGED) &&
((fabs(pcoords[1]-params[1])) < CONVERGED) &&
((fabs(pcoords[2]-params[2])) < CONVERGED) )
{
converged = 1;
}
//
// if not converged, repeat
//
else
{
params[0] = pcoords[0];
params[1] = pcoords[1];
params[2] = pcoords[2];
}
}
//
// if not converged, set the parametric coordinates to arbitrary values
// outside of element
//
if ( !converged )
{
pcoords[0] = pcoords[1] = pcoords[2] = 10.0;
return LARGE_FLOAT;
}
else
{
if ( pcoords[0] >= -1.0 && pcoords[1] <= 1.0 &&
pcoords[1] >= -1.0 && pcoords[1] <= 1.0 &&
pcoords[2] >= -1.0 && pcoords[2] <= 1.0 )
{
for(i=0; i<3; i++) pcoords[i] = 0.5*(pcoords[i]+1.0); // shift to (0,1)
return 0.0; // inside hexahedron
}
else
{
for (i=0; i<3; i++)
{
if (pcoords[i] < -1.0) pcoords[i] = -1.0;
if (pcoords[i] > 1.0) pcoords[i] = 1.0;
}
this->EvaluateLocation(subId, pcoords, closestPoint);
for(i=0; i<3; i++) pcoords[i] = 0.5*(pcoords[i]+1.0); // shift to (0,1)
return math.Distance2BetweenPoints(closestPoint,x);
}
}
}
//
// Compute iso-parametrix shape functions
//
void vlHexahedron::ShapeFunctions(float pcoords[3], float sf[8])
{
double rm, rp, sm, sp, tm, tp;
rm = 1. - pcoords[0];
rp = 1. + pcoords[0];
sm = 1. - pcoords[1];
sp = 1. + pcoords[1];
tm = 1. - pcoords[2];
tp = 1. + pcoords[2];
sf[0] = 0.125*rm*sm*tm;
sf[1] = 0.125*rp*sm*tm;
sf[2] = 0.125*rp*sp*tm;
sf[3] = 0.125*rm*sp*tm;
sf[4] = 0.125*rm*sm*tp;
sf[5] = 0.125*rp*sm*tp;
sf[6] = 0.125*rp*sp*tp;
sf[7] = 0.125*rm*sp*tp;
}
void vlHexahedron::ShapeDerivs(float pcoords[3], float derivs[24])
{
double rm, rp, sm, sp, tm, tp;
rm = 1. - pcoords[0];
rp = 1. + pcoords[0];
sm = 1. - pcoords[1];
sp = 1. + pcoords[1];
tm = 1. - pcoords[2];
tp = 1. + pcoords [2];
derivs[0] = -0.125*sm*tm;
derivs[1] = 0.125*sm*tm;
derivs[2] = 0.125*sp*tm;
derivs[3] = -0.125*sp*tm;
derivs[4] = -0.125*sm*tp;
derivs[5] = 0.125*sm*tp;
derivs[6] = 0.125*sp*tp;
derivs[7] = -0.125*sp*tp;
derivs[8] = -0.125*rm*tm;
derivs[9] = -0.125*rp*tm;
derivs[10] = 0.125*rp*tm;
derivs[11] = 0.125*rm*tm;
derivs[12] = -0.125*rm*tp;
derivs[13] = -0.125*rp*tp;
derivs[14] = 0.125*rp*tp;
derivs[15] = 0.125*rm*tp;
derivs[16] = -0.125*rm*sm;
derivs[17] = -0.125*rp*sm;
derivs[18] = -0.125*rp*sp;
derivs[19] = -0.125*rm*sp;
derivs[20] = 0.125*rm*sm;
derivs[21] = 0.125*rp*sm;
derivs[22] = 0.125*rp*sp;
derivs[23] = 0.125*rm*sp;
}
void vlHexahedron::EvaluateLocation(int& subId, float pcoords[3], float x[3])
{
int i, j;
float sf[8], *pt, pc[3];
for (i=0;i<3;i++) pc[i] = 2.0*pcoords[i] - 1.0; //shift to -1<=r,s,t<=1
this->ShapeFunctions(pc, sf);
x[0] = x[1] = x[2] = 0.0;
for (i=0; i<8; i++)
{
pt = this->Points.GetPoint(i);
for (j=0; j<3; j++)
{
x[j] += pt[j] * sf[i];
}
}
}
/*=========================================================================
Program: Visualization Library
Module: Pixel.cc
Language: C++
Date: $Date$
Version: $Revision$
Description:
---------------------------------------------------------------------------
This file is part of the Visualization Library. No part of this file
or its contents may be copied, reproduced or altered in any way
without the express written consent of the authors.
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen 1993, 1994
=========================================================================*/
#include <math.h>
#include "Rect.hh"
#include "Polygon.hh"
#include "Plane.hh"
#include "vlMath.hh"
//
// Note: the ordering of the Points and PointIds is important. See text.
//
float vlRectangle::EvaluatePosition(float x[3], int& subId, float pcoords[3])
{
float *pt1, *pt2, *pt3;
vlPolygon poly;
vlPlane plane;
int i;
float xProj[3], closestPoint[3], p[3], p21[3], p31[3];
float l21, l31, n[3];
vlMath math;
subId = 0;
pcoords[0] = pcoords[1] = pcoords[2] = 0.0;
//
// Get normal for rectangle
//
pt1 = this->Points.GetPoint(0);
pt2 = this->Points.GetPoint(1);
pt3 = this->Points.GetPoint(2);
poly.ComputeNormal (pt1, pt2, pt3, n);
//
// Project point to plane
//
plane.ProjectPoint(x,pt1,n,xProj);
for (i=0; i<3; i++)
{
p21[i] = pt2[i] - pt1[i];
p31[i] = pt3[i] - pt1[i];
p[i] = x[i] - pt1[i];
}
if ( (l21=math.Norm(p21)) == 0.0 ) l21 = 1.0;
if ( (l31=math.Norm(p31)) == 0.0 ) l31 = 1.0;
pcoords[0] = math.Dot(p21,p) / l21;
pcoords[1] = math.Dot(p31,p) / l31;
if ( pcoords[0] >= -1.0 && pcoords[1] <= 1.0 &&
pcoords[1] >= -1.0 && pcoords[1] <= 1.0 )
{
return math.Distance2BetweenPoints(xProj,x); //projection distance
}
else
{
for (i=0; i<2; i++)
{
if (pcoords[i] < -1.0) pcoords[i] = -1.0;
if (pcoords[i] > 1.0) pcoords[i] = 1.0;
}
this->EvaluateLocation(subId, pcoords, closestPoint);
return math.Distance2BetweenPoints(closestPoint,x);
}
}
void vlRectangle::EvaluateLocation(int& subId, float pcoords[3], float x[3])
{
float *pt1, *pt2, *pt3;
int i;
pt1 = this->Points.GetPoint(0);
pt2 = this->Points.GetPoint(1);
pt3 = this->Points.GetPoint(2);
for (i=0; i<3; i++)
{
x[i] = pt1[i] + pcoords[0]*(pt2[i] - pt1[i]) +
pcoords[1]*(pt3[i] - pt1[i]);
}
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment