Commit f5f0c3da authored by allens's avatar allens
Browse files

moved coordinate conversion to the fields and inital implementation for NIMROD interpolation

git-svn-id: http://visit.ilight.com/svn/visit/trunk/src@14788 18c085ea-50e0-402c-830e-de6fd14e8384
parent f6dd1395
......@@ -56,7 +56,6 @@ avtIVPSolver.C
avtIVPM3DC1Field.C
avtIVPM3DC1Integrator.C
avtIVPNIMRODField.C
avtIVPNIMRODIntegrator.C
avtIVPVTKField.C
avtIVPVTKTimeVaryingField.C
avtPoincareIC.C
......
......@@ -95,7 +95,6 @@ avtIVPAdamsBashforth::avtIVPAdamsBashforth()
initialized = 0;
degenerate_iterations = 0;
stiffness_eps = tol / 1000.0;
convertToCartesian = false;
}
// ****************************************************************************
......@@ -507,8 +506,8 @@ avtIVPAdamsBashforth::Step(avtIVPField* field, double t_max,
if( convertToCartesian )
{
(*ivpstep)[0] = CylindricalToCartesian( yCur );
(*ivpstep)[1] = CylindricalToCartesian( yNew );
(*ivpstep)[0] = field->ConvertToCartesian( yCur );
(*ivpstep)[1] = field->ConvertToCartesian( yNew );
}
else
{
......
......@@ -104,8 +104,6 @@ avtIVPDopri5::avtIVPDopri5()
h_max = 0.0;
nonsti = 0;
convertToCartesian = false;
}
......@@ -689,15 +687,15 @@ avtIVPDopri5::Step(avtIVPField* field, double t_max,
if( convertToCartesian )
{
(*ivpstep)[0] = CylindricalToCartesian( y );
(*ivpstep)[1] = CylindricalToCartesian( y + (h*k1/4.) );
(*ivpstep)[2] = CylindricalToCartesian( (y + y_new)/2 +
(*ivpstep)[0] = field->ConvertToCartesian( y );
(*ivpstep)[1] = field->ConvertToCartesian( y + (h*k1/4.) );
(*ivpstep)[2] = field->ConvertToCartesian( (y + y_new)/2 +
h*( (d1+1)*k1 +
d3*k3 + d4*k4 +
d5*k5 + d6*k6 +
(d7-1)*k7 )/6. );
(*ivpstep)[3] = CylindricalToCartesian( y_new - h*k7/4 );
(*ivpstep)[4] = CylindricalToCartesian( y_new );
(*ivpstep)[3] = field->ConvertToCartesian( y_new - h*k7/4 );
(*ivpstep)[4] = field->ConvertToCartesian( y_new );
}
else
{
......@@ -783,5 +781,3 @@ avtIVPDopri5::AcceptStateVisitor(avtIVPStateHelper& aiss)
.Accept(y)
.Accept(k1);
}
......@@ -103,6 +103,10 @@ class IVP_API avtIVPField
virtual avtVector operator()(const double& t,
const avtVector& x) const = 0;
virtual avtVector ConvertToCartesian(const avtVector& pt) const = 0;
virtual avtVector ConvertToCylindrical(const avtVector& pt) const = 0;
virtual double ComputeVorticity(const double& t,
const avtVector& x ) const = 0;
virtual double ComputeScalarVariable(unsigned char index,
......
......@@ -782,6 +782,41 @@ avtIVPM3DC1Field::operator()( const double &t, const avtVector &p ) const
}
// ****************************************************************************
// Method: avtIVPM3DC1Field::ConvertToCartesian
//
// Purpose: Converts the coordinates from local cylindrical to
// cartesian coordinates
//
// Programmer: Christoph Garth
// Creation: February 25, 2008
//
// ****************************************************************************
avtVector
avtIVPM3DC1Field::ConvertToCartesian(const avtVector& pt) const
{
return avtVector(pt[0]*cos(pt[1]), pt[0]*sin(pt[1]), pt[2] );
}
// ****************************************************************************
// Method: avtIVPM3DC1Field::ConvertToCylindrical
//
// Purpose: Converts the coordinates from local cylindrical to
// cartesian coordinates
//
// Programmer: Christoph Garth
// Creation: February 25, 2008
//
// ****************************************************************************
avtVector
avtIVPM3DC1Field::ConvertToCylindrical(const avtVector& pt) const
{
return pt;
}
// ****************************************************************************
// Method: interpBcomps
//
......
......@@ -99,6 +99,9 @@ class IVP_API avtIVPM3DC1Field: public avtIVPVTKField
avtVector operator()( const double &t, const avtVector &v ) const;
avtVector ConvertToCartesian(const avtVector& pt) const;
avtVector ConvertToCylindrical(const avtVector& pt) const;
void interpBcomps(float *B, double *x, int element, double *xieta) const;
float interp(float *var, int el, double *lcoords) const;
......
......@@ -44,85 +44,58 @@
#include <DebugStream.h>
#include <vtkCellData.h>
#include <vtkStructuredGrid.h>
#include <vtkIntArray.h>
#include <vtkFloatArray.h>
#include <vtkIdList.h>
#include <math.h>
#define ELEMENT_SIZE 7
#define SCALAR_SIZE 20
// ****************************************************************************
// Method: avtIVPNIMRODField constructor
//
// Creationist: Allen Sanderson
// Creation: 20 November 2009
//
// Modifications:
//
// Hank Childs, Mon Jun 7 14:48:03 CDT 2010
// Initialize variables to prevent compiler warning. (They are unused and
// just to get the template instantiation right.)
//
// ****************************************************************************
avtIVPNIMRODField::avtIVPNIMRODField( vtkDataSet* dataset,
avtCellLocator* locator ) :
avtIVPVTKField( dataset, locator )
avtCellLocator* locator ) :
avtIVPVTKField( dataset, locator ),
grid_fourier_series( 0 ), data_fourier_series( 0 ),
Drad( 2 ), Dtheta( 2 )
{
vtkStructuredGrid* grid = vtkStructuredGrid::SafeDownCast( dataset );
// Pick off all of the data stored with the vtk field.
// Get the numver of elements for checking the validity of the data.
if( grid->GetDataDimension() == 3 )
{
int dims[3];
// Because the triangluar mesh is defined by using non unique points
// and the data is cell centered data VisIt moves it out to the
// nodes for STREAMLINES thus there are 3 times the number of
// original values.
if( ds->GetPointData()->GetArray("hidden/elements") )
nelms =
ds->GetPointData()->GetArray("hidden/elements")->GetNumberOfTuples() / 3;
grid->GetDimensions( dims );
// 2.0 Change data is at the cells for POINCARE
Nrad = dims[0];
Ntheta = dims[1];
Nphi = dims[2];
}
else
nelms =
ds->GetCellData()->GetArray("hidden/elements")->GetNumberOfTuples();
{
debug1 << "Mesh dimension is not 3 " << endl;
return;
}
// Dummy variable to the template class
int *intPtr, intVar = 0;
float *fltPtr, fltVar = 0;
// Single values from the header attributes.
intPtr = SetDataPointer( ds, intVar, "hidden/header/linear", nelms, 1 );
linflag = intPtr[0];
delete [] intPtr;
intPtr = SetDataPointer( ds, intVar, "hidden/header/ntor", nelms, 1 );
tmode = intPtr[0];
delete [] intPtr;
fltPtr = SetDataPointer( ds, fltVar, "hidden/header/bzero", nelms, 1 );
bzero = fltPtr[0];
delete [] fltPtr;
fltPtr = SetDataPointer( ds, fltVar, "hidden/header/rzero", nelms, 1 );
rzero = fltPtr[0];
delete [] fltPtr;
float fltVar = 0;
// The mesh elements.
elements = SetDataPointer( ds, fltVar, "hidden/elements", nelms, ELEMENT_SIZE );
// grid_fourier_series =
// SetDataPointer( ds, fltVar, "hidden/grid_fourier_series", 3 );
// Vector values from the field.
psi0 = SetDataPointer( ds, fltVar, "hidden/equilibrium/psi", nelms, SCALAR_SIZE );
f0 = SetDataPointer( ds, fltVar, "hidden/equilibrium/f", nelms, SCALAR_SIZE );
psinr = SetDataPointer( ds, fltVar, "hidden/psi", nelms, SCALAR_SIZE );
psini = SetDataPointer( ds, fltVar, "hidden/psi_i", nelms, SCALAR_SIZE );
fnr = SetDataPointer( ds, fltVar, "hidden/f", nelms, SCALAR_SIZE );
fni = SetDataPointer( ds, fltVar, "hidden/f_i", nelms, SCALAR_SIZE );
// Now set some values using the above data.
F0 = -bzero * rzero;
findElementNeighbors();
// data_fourier_series =
// SetDataPointer( ds, fltVar, "hidden/data_fourier_series", 3 );
}
// ****************************************************************************
......@@ -133,11 +106,30 @@ avtIVPNIMRODField::avtIVPNIMRODField( vtkDataSet* dataset,
//
// ****************************************************************************
avtIVPNIMRODField::avtIVPNIMRODField( float *elementsPtr, int nelements )
: avtIVPVTKField(NULL, NULL), elements( elementsPtr), neighbors(0),
psi0(0), f0(0), psinr(0), psini(0), fnr(0), fni(0), nelms(nelements)
avtIVPNIMRODField::avtIVPNIMRODField( unsigned int nRad,
unsigned int nTheta,
unsigned int nPhi,
float *gfs,
float *dfs )
: avtIVPVTKField(NULL, NULL),
grid_fourier_series( gfs ), data_fourier_series( dfs ),
Nrad( nRad ), Ntheta( nTheta ), Nphi( nPhi ),
Drad( 2 ), Dtheta( 2 )
{
findElementNeighbors();
lagrange_nodes[0][0] = 0.0;
lagrange_nodes[1][0] = 0.0; lagrange_nodes[1][1] = 1.0;
lagrange_nodes[2][0] = 0.0; lagrange_nodes[2][1] = 0.5; lagrange_nodes[2][2] = 0.5;
lagrange_nodes[3][0] = 0.0; lagrange_nodes[3][1] = 0.276393202250021031; lagrange_nodes[3][2] = 0.723606797749978969; lagrange_nodes[3][3] = 1.0;
lagrange_nodes[4][0] = 0.0; lagrange_nodes[4][1] = 0.172673164646011429; lagrange_nodes[4][2] = 0.500000000000000000; lagrange_nodes[4][3] = 0.827326835353988571; lagrange_nodes[4][4] = 1.0;
lagrange_nodes[5][0] = 0.0; lagrange_nodes[5][1] = 0.117472338035267654; lagrange_nodes[5][2] = 0.357384241759677453; lagrange_nodes[5][3] = 0.642615758240322547; lagrange_nodes[5][4] = 0.882527661964732346; lagrange_nodes[5][5] = 1.0;
// { 0.0 },
// { 0.0, 1.0 },
// { 0.0, 0.5, 1.0 },
// { 0.0, 0.276393202250021031, 0.723606797749978969, 1.0 },
// { 0.0, 0.172673164646011429, 0.500000000000000000, 0.827326835353988571, 1.0 },
// { 0.0, 0.117472338035267654, 0.357384241759677453, 0.642615758240322547, 0.882527661964732346, 1.0 },
// };
}
......@@ -151,16 +143,8 @@ avtIVPNIMRODField::avtIVPNIMRODField( float *elementsPtr, int nelements )
avtIVPNIMRODField::~avtIVPNIMRODField()
{
if( neighbors ) free(neighbors);
if( trigtable ) free(trigtable);
if( elements ) delete [] elements;
if( psi0 ) delete [] psi0;
if( f0 ) delete [] f0;
if( psinr ) delete [] psinr;
if( psini ) delete [] psini;
if( fnr ) delete [] fnr;
if( fni ) delete [] fni;
if( grid_fourier_series ) free(grid_fourier_series);
if( data_fourier_series ) free(data_fourier_series);
}
......@@ -176,7 +160,6 @@ template< class type >
type* avtIVPNIMRODField::SetDataPointer( vtkDataSet *ds,
const type var,
const char* varname,
const int ntuples,
const int ncomponents )
{
vtkDataArray *array;
......@@ -189,13 +172,6 @@ type* avtIVPNIMRODField::SetDataPointer( vtkDataSet *ds,
if( ds->GetPointData()->GetArray(varname) )
{
array = ds->GetPointData()->GetArray(varname);
XX = 3;
}
// 2.0 Change data is at the cells for POINCARE
else
{
array = ds->GetCellData()->GetArray(varname);
XX = 1;
}
if( array == 0 )
......@@ -207,7 +183,9 @@ type* avtIVPNIMRODField::SetDataPointer( vtkDataSet *ds,
return 0;
}
if( ntuples != array->GetNumberOfTuples() / XX ||
int ntuples = Nrad*Ntheta*Nphi;
if( ntuples != array->GetNumberOfTuples() ||
ncomponents != array->GetNumberOfComponents() )
{
debug1 << "Variable " << varname
......@@ -217,12 +195,6 @@ type* avtIVPNIMRODField::SetDataPointer( vtkDataSet *ds,
return 0;
}
// 2.0 Change data is no longer at the points but still is at the
// cells so we should be able to use the pointer directly but for
// some reason it causes problems.
//
// return (type*) array->GetVoidPointer(0);
type* newptr = new type[ntuples*ncomponents];
if( newptr == 0 )
......@@ -231,19 +203,13 @@ type* avtIVPNIMRODField::SetDataPointer( vtkDataSet *ds,
return 0;
}
// Because the triangluar mesh is defined by using non unique points
// and the data is cell centered data VisIt moves it out to the
// nodes. So create a new structure that is what is really needed.
// 2.0 Change data is no longer at the points but still is at the
// cells so the above is no longer valid.
if( array->IsA("vtkIntArray") )
{
int* ptr = (int*) array->GetVoidPointer(0);
for( int i=0; i<ntuples; ++i )
for( int j=0; j<ncomponents; ++j )
newptr[i*ncomponents+j] = ptr[i*XX*ncomponents+j];
newptr[i*ncomponents+j] = ptr[i*ncomponents+j];
return newptr;
}
......@@ -253,7 +219,7 @@ type* avtIVPNIMRODField::SetDataPointer( vtkDataSet *ds,
for( int i=0; i<ntuples; ++i )
for( int j=0; j<ncomponents; ++j )
newptr[i*ncomponents+j] = ptr[i*XX*ncomponents+j];
newptr[i*ncomponents+j] = ptr[i*ncomponents+j];
return newptr;
}
......@@ -263,7 +229,7 @@ type* avtIVPNIMRODField::SetDataPointer( vtkDataSet *ds,
for( int i=0; i<ntuples; ++i )
for( int j=0; j<ncomponents; ++j )
newptr[i*ncomponents+j] = ptr[i*XX*ncomponents+j];
newptr[i*ncomponents+j] = ptr[i*ncomponents+j];
return newptr;
}
......@@ -286,9 +252,6 @@ type* avtIVPNIMRODField::SetDataPointer( vtkDataSet *ds,
// Gets the B field components directly - should not be used for
// calculating integral curves.
//
// THIS CODE SHOULD NOT BE USED FOR FIELDLINE INTEGRATION!!!!
//
//
// Programmer: Allen Sanderson
// Creation: October 24, 2009
//
......@@ -299,549 +262,284 @@ type* avtIVPNIMRODField::SetDataPointer( vtkDataSet *ds,
avtVector
avtIVPNIMRODField::operator()( const double &t, const avtVector &p ) const
{
// NOTE: Assumes the point is in cylindrical coordiantes.
double pt[3] = { p[0], p[1], p[2] };
avtVector linear_vec = avtIVPVTKField::operator()( t, p );
pt[0] = p[0];
pt[1] = p[1];
pt[2] = p[2];
if( !FindCell( t, p ) )
throw Undefined();
/* Find the element containing the point; get local coords xi,eta */
double xieta[2];
int element;
vtkIdList *ptIds = vtkIdList::New();
if ((element = get_tri_coords2D(pt, xieta)) < 0)
{
return avtVector(0,0,0);
}
else
{
float B[3];
ds->GetCellPoints( lastCell, ptIds );
interpBcomps(B, pt, element, xieta);
double center[3], pt[3];
// The B value is in cylindrical coordiantes
return avtVector( B[0], B[1], B[2] );
for( unsigned int i=0; i<Nphi; ++ i )
{
ds->GetPoint( i+62*101, center);
cerr << i << " "
<< center[0] << " " << center[1] << " " << center[2] << " "
<< sqrt(center[0]*center[0]+center[1]*center[1]) << " " << endl;
}
// converstion to a right hand system.
// avtVector( B[0] * cos(pt[1]) - B[1] * sin(pt[1]),
// B[0] * sin(pt[1]) + B[1] * cos(pt[1]),
// B[2] );
}
// ****************************************************************************
// Method: avtIVPNIMRODField IsInside
//
// Creationist: Allen Sanderson
// Creation: 16 April 2010
//
// The VTK check should work but is not reliable because the mesh
// lies in the XZ plane and Y may or may not be exactly zero. As such
// due to round off VTK may say the point is outside the
// mesh. As such, use the native check instead.
//
// ****************************************************************************
cerr << lastCell << " ( " << endl;
bool avtIVPNIMRODField::IsInside(const double& t, const avtVector& x) const
{
double xin[3], xout[3];
unsigned int iRad=Nrad, iTheta=Ntheta, iPhi=Nphi;
xin[0] = x[0];
xin[1] = x[1];
xin[2] = x[2];
double minRad=+1e10, minTheta=+1e10, minPhi=+1e10;
double maxRad=-1e10, maxTheta=-1e10, maxPhi=-1e10;
return (bool) get_tri_coords2D(xin, xout);
}
unsigned int ic=0;
// ****************************************************************************
// Method: findElementNeighbors
//
// Creationist: Allen Sanderson
// Creation: 20 November 2009
//
// ****************************************************************************
void avtIVPNIMRODField::findElementNeighbors()
{
v_entry *vert_list = 0;
edge *edge_list = 0;
float *ptr;
double x[3], y[3], co, sn;
int el, vert, tri[3], vlen;
/* Allocate, initialize neighbor table */
neighbors = (int *)malloc(3 * nelms * sizeof(int));
if (neighbors == NULL) {
fputs("Insufficient memory in findElementNeighbors.\n", stderr);
exit(1);
}
for( unsigned int i=0; i<ptIds->GetNumberOfIds(); ++i )
{
unsigned int tRad, tTheta, tPhi, tmp = ptIds->GetId(i);
for (el=0; el<3*nelms; el++)
neighbors[el] = -1;
cerr << tmp << " ";
/* Allocate trig table */
trigtable = (double *)malloc(2 * nelms * sizeof(double));
if (trigtable == NULL) {
fputs("Insufficient memory in findElementNeighbors.\n", stderr);
exit(1);
}
tRad = tmp / (Ntheta*Nphi);
/* Allocate, initialize temporary hash tables */
vert_list = (v_entry *)malloc(3 * nelms * sizeof(v_entry));
if (vert_list == NULL) {
fputs("Insufficient memory in findElementNeighbors.\n", stderr);
exit(1);
}
vlen = 0;
edge_list = (edge *)malloc(3 * nelms * sizeof(edge));
if (edge_list == NULL) {
fputs("Insufficient memory in findElementNeighbors.\n", stderr);
exit(1);
}
tmp = tmp % (Ntheta*Nphi);
for (vert=0; vert<3*nelms; vert++)
edge_list[vert].n = 0;
tTheta = tmp / Nphi;
/* Loop over elements, finding vertices, edges, neighbors */
tPhi = tmp % Nrad;
if( iRad >= tRad )
{
iRad = tRad;
//For each element, the first 6 values are a, b, c, theta, x, and z.
//The nodes of the element are located at
// (x,z),
// (x+(a+b)*cos(theta),z+(a+b)*sin(theta)),
// (x+b*cos(theta)-c*sin(theta),z+b*sin(theta)+c*cos(theta)).
if( iTheta >= tTheta )
{
iTheta = tTheta;
for (el=0; el<nelms; el++) {
ptr = elements + ELEMENT_SIZE*el;
co = trigtable[2*el] = cos(ptr[3]);
sn = trigtable[2*el + 1] = sin(ptr[3]);
if( iPhi > tPhi )
{
iPhi = tPhi;
ic = i;
}
}
}
x[0] = ptr[4];
y[0] = ptr[5];
x[1] = x[0] + (ptr[0] + ptr[1])*co;
y[1] = y[0] + (ptr[0] + ptr[1])*sn;
ds->GetPoint( ptIds->GetId(i), pt);
ds->GetPoint( tPhi, center);
x[2] = x[0] + ptr[1]*co - ptr[2]*sn;
y[2] = y[0] + ptr[1]*sn + ptr[2]*co;
cerr << " " << tPhi << " " << tTheta << " " << tRad << " "
// << center[0] << " " << center[1] << " " << center[2] << " ";
<< pt[0] << " " << pt[1] << " " << pt[2] << " ";
for (vert=0; vert<3; vert++)
register_vert(vert_list, &vlen, x[vert], y[vert], tri+vert);
double dx = pt[0] - center[0];
double dz = pt[2] - center[2];
for (vert=0; vert<3; vert++)
add_edge(edge_list, tri, vert, el, neighbors);
} /* end loop el */
double rad = sqrt( dx*dx + dz*dz );
double theta = atan2( dz, dx );
double phi = atan2( pt[1], pt[0] );
// fprintf(stderr, "%d / %d unique vertices\n", vlen, 3*nelms);
// fprintf(stderr, "Neighbors of element 0: %d, %d, %d\n", neighbors[0],
// neighbors[1], neighbors[2]);
cerr << phi << " " << theta << " " << rad << " " << endl;
/* Use unique vert list to find mesh bounds */
Rmin = Rmax = vert_list[0].x;
zmin = zmax = vert_list[0].y;
for (vert=1; vert<vlen; vert++) {
if (Rmin > vert_list[vert].x) Rmin = vert_list[vert].x;
if (Rmax < vert_list[vert].x) Rmax = vert_list[vert].x;
if (zmin > vert_list[vert].y) zmin = vert_list[vert].y;
if (zmax < vert_list[vert].y) zmax = vert_list[vert].y;
}
if( minRad > rad ) minRad = rad;
if( maxRad < rad ) maxRad = rad;
// fprintf(stderr, "R bounds: %lf, %lf\nz bounds: %lf, %lf\n",
// Rmin, Rmax, zmin, zmax);
free(vert_list);
free(edge_list);
}
if( minTheta > theta ) minTheta = theta;
if( maxTheta < theta ) maxTheta = theta;
if( minPhi > phi ) minPhi = phi;
if( maxPhi < phi ) maxPhi = phi;