Commit 45053410 authored by allens's avatar allens

added new class for fieldline

git-svn-id: http://visit.ilight.com/svn/visit/trunk/src@16708 18c085ea-50e0-402c-830e-de6fd14e8384
parent add50505
......@@ -111,7 +111,8 @@ enum FieldlineType { UNKNOWN_TYPE = 0,
IRRATIONAL = 20,
FLUX_SURFACE = 21,
ISLAND_CHAIN = 22,
ISLANDS_WITHIN_ISLANDS = 23,
ISLAND_WITH_SECONDARY_ISLANDS = 23,
ISLAND_AMBIGUOUS_AXIS = 24,
CHAOTIC = 30 };
......
......@@ -449,13 +449,14 @@ FieldlineLib::convexHull( std::vector< std::pair< Point, unsigned int > > &hullP
// Stop when the first point is found again.
} while( min != npts );
// for( unsigned int i=0; i<m; ++i ) {
// if( verboseFlag )
// std::cerr << "convexHull " << hullPts[i].second << std::endl;
// }
for( unsigned int i=0; i<m; ++i ) {
if( verboseFlag )
std::cerr << "convexHull " << hullPts[i].second << " "
<< hullPts[i].first << std::endl;
}
// if( verboseFlag )
// std::cerr << std::endl;
if( verboseFlag )
std::cerr << std::endl;
}
......@@ -916,7 +917,7 @@ unsigned int FieldlineLib::Blankinship( unsigned int toroidalWinding,
unsigned int offset = GCD( toroidalWinding, poloidalWinding );
if( toroidalWinding == poloidalWinding )
skip = 1;
skip = 0;
else if( toroidalWinding > 1 && poloidalWinding != 0 ) {
// To find the skip find the mutual primes via the
......@@ -2320,6 +2321,10 @@ FieldlineLib::fieldlineProperties( std::vector< Point > &ptList,
}
bool helicity = ccw( ptList[1] - ptList[0], ptList[1] - ptList[2] );
std::cerr << "helicity " << helicity << std::endl;
// Get the safety factor.
// for( unsigned int i=0; i<poloidal_puncture_pts.size(); ++i )
// std::cerr << i << " "
......@@ -2849,7 +2854,7 @@ FieldlineLib::fieldlineProperties( std::vector< Point > &ptList,
// the toroidal resonance
if( (type == FieldlineProperties::UNKNOWN_TYPE ||
type == FieldlineProperties::ISLAND_CHAIN ||
type == FieldlineProperties::ISLANDS_WITHIN_ISLANDS ) &&
type == FieldlineProperties::ISLAND_WITH_SECONDARY_ISLANDS ) &&
( (toroidalResonance > 1 /* && poloidalResonance >= 1 */) || // Always true.
......@@ -2861,8 +2866,8 @@ FieldlineLib::fieldlineProperties( std::vector< Point > &ptList,
// of islands.
if( toroidalWinding == poloidalWinding )
{
toroidalWinding = 1;
poloidalWinding = 1;
// toroidalWinding = 1;
// poloidalWinding = 1;
poloidalWindingP = poloidalWinding;
}
else
......@@ -2892,7 +2897,7 @@ FieldlineLib::fieldlineProperties( std::vector< Point > &ptList,
// smaller islands around an island.
else
{
type = FieldlineProperties::ISLANDS_WITHIN_ISLANDS;
type = FieldlineProperties::ISLAND_WITH_SECONDARY_ISLANDS;
if( verboseFlag )
std::cerr << "Secondary resonances = "
......@@ -2927,7 +2932,7 @@ FieldlineLib::fieldlineProperties( std::vector< Point > &ptList,
unsigned int nodes = poloidal_puncture_pts.size() / islands / 2;
// Add two more puncture points per island within and island.
if( type == FieldlineProperties::ISLANDS_WITHIN_ISLANDS )
if( type == FieldlineProperties::ISLAND_WITH_SECONDARY_ISLANDS )
nPuncturesNeeded = (nodes + 2) * islands * 2;
// Add five more puncture points per island.
else //if( type == FieldlineProperties::ISLAND_CHAIN )
......@@ -3103,7 +3108,7 @@ FieldlineLib::fieldlineProperties( std::vector< Point > &ptList,
unsigned int m = 0; // starting index
convexHull( hullPts, m, toroidalWinding, 1 );
convexHull( hullPts, m, toroidalWinding, helicity );
if( m != toroidalWinding )
{
......@@ -3121,9 +3126,11 @@ FieldlineLib::fieldlineProperties( std::vector< Point > &ptList,
// Find all the differences and count each one.
for(unsigned int i=0; i<toroidalWinding; ++i )
{
offset = ( (hullPts[i ].second -
hullPts[(i+1)%toroidalWinding].second) +
toroidalWinding ) % toroidalWinding;
unsigned int i1 = (i+1) % toroidalWinding;
// Offset from one puncture point to it's neighbor.
offset = (hullPts[i1].second - hullPts[i].second + toroidalWinding) %
toroidalWinding;
// Find this offset in the list.
ic = offsets.find( offset );
......@@ -3131,7 +3138,7 @@ FieldlineLib::fieldlineProperties( std::vector< Point > &ptList,
// Not found, new offset.
if( ic == offsets.end() )
offsets.insert( std::pair<int, int>( offset, 1) );
// Found this difference, increment the count.
// Found this offset, increment the count.
else
(*ic).second++;
}
......@@ -3140,7 +3147,7 @@ FieldlineLib::fieldlineProperties( std::vector< Point > &ptList,
if( verboseFlag )
std::cerr << "Multiple offsets ";
// Find the difference that occurs most often.
// Find the offset that occurs most often.
unsigned int nMatches = 0;
ic = offsets.begin();
......@@ -3164,10 +3171,16 @@ FieldlineLib::fieldlineProperties( std::vector< Point > &ptList,
if( verboseFlag )
std::cerr << std::endl;
// Secondary angle around the nonaxisymmetric island.
// Secondary rotation around the nonaxisymmetric island.
if( offset != 1 && offset != toroidalWinding-1 )
{
poloidalWindingP = Blankinship( toroidalWinding, offset );
windingGroupOffset = offset;
nnodes = 1;
type = FieldlineProperties::ISLAND_AMBIGUOUS_AXIS;
if( verboseFlag )
std::cerr << "Secondary poloidal rotation "
......@@ -3425,7 +3438,8 @@ FieldlineLib::fieldlineProperties( std::vector< Point > &ptList,
double local_safetyFactor =
(double) toroidalWinding / (double) poloidalWinding;
windingGroupOffset = Blankinship( toroidalWinding, poloidalWinding );
if( poloidalWinding == poloidalWindingP )
windingGroupOffset = Blankinship( toroidalWinding, poloidalWinding );
nnodes = poloidal_puncture_pts.size() / toroidalWinding;
......@@ -3634,7 +3648,7 @@ FieldlineLib::fieldlineProperties2( std::vector< Point > &ptList,
// only check the toroidalFirstResonance
if( (type == FieldlineProperties::UNKNOWN_TYPE ||
type == FieldlineProperties::ISLAND_CHAIN ||
type == FieldlineProperties::ISLANDS_WITHIN_ISLANDS ) &&
type == FieldlineProperties::ISLAND_WITH_SECONDARY_ISLANDS ) &&
toroidalResonance > 1 ) // && poloidalResonance >= 1 ) // Always true.
{
......@@ -3661,7 +3675,7 @@ FieldlineLib::fieldlineProperties2( std::vector< Point > &ptList,
// smaller islands around an island.
else
{
type = FieldlineProperties::ISLANDS_WITHIN_ISLANDS;
type = FieldlineProperties::ISLAND_WITH_SECONDARY_ISLANDS;
if( verboseFlag )
std::cerr << "Secondary resonances = "
......@@ -3692,7 +3706,7 @@ FieldlineLib::fieldlineProperties2( std::vector< Point > &ptList,
unsigned int nodes = poloidal_puncture_pts.size() / islands / 2;
// Add two more puncture points per island within and island.
if( type == FieldlineProperties::ISLANDS_WITHIN_ISLANDS )
if( type == FieldlineProperties::ISLAND_WITH_SECONDARY_ISLANDS )
nPuncturesNeeded = (nodes + 2) * islands * 2;
// Add five more puncture points per island.
else //if( type == FieldlineProperties::ISLAND_CHAIN )
......
......@@ -112,7 +112,8 @@ enum FieldlineType { UNKNOWN_TYPE = 0,
IRRATIONAL = 20,
FLUX_SURFACE = 21,
ISLAND_CHAIN = 22,
ISLANDS_WITHIN_ISLANDS = 23,
ISLAND_WITH_SECONDARY_ISLANDS = 23,
ISLAND_AMBIGUOUS_AXIS = 24,
CHAOTIC = 30 };
......
......@@ -405,7 +405,10 @@ avtPoincareFilter::ContinueExecute()
FieldlineProperties::ISLAND_CHAIN ||
poincare_ic->properties.type ==
FieldlineProperties::ISLANDS_WITHIN_ISLANDS) &&
FieldlineProperties::ISLAND_WITH_SECONDARY_ISLANDS ||
poincare_ic->properties.type ==
FieldlineProperties::ISLAND_AMBIGUOUS_AXIS) &&
poincare_ic->properties.analysisState ==
FieldlineProperties::ADD_O_POINTS &&
......@@ -768,7 +771,7 @@ avtPoincareFilter::CreatePoincareOutput( avtDataTree *dt,
unsigned int islandGroups = properties.islandGroups;
unsigned int nnodes = properties.nnodes;
std::vector< avtVector > &OPoints = properties.OPoints;
std::vector< avtVector > &OPoints = properties.OPoints;
bool completeIslands = true;
......@@ -788,11 +791,6 @@ avtPoincareFilter::CreatePoincareOutput( avtDataTree *dt,
<< toroidalWinding << "," << poloidalWinding << " ("
<< safetyFactor << ") ";
if( poloidalWinding != poloidalWindingP )
std::cerr << toroidalWinding << "," << poloidalWindingP << " ("
<< (float) toroidalWinding / (float) poloidalWindingP << ") ";
if( type == FieldlineProperties::RATIONAL )
std::cerr << "rational surface ";
......@@ -803,11 +801,19 @@ avtPoincareFilter::CreatePoincareOutput( avtDataTree *dt,
std::cerr << islands << " island chain with resonances: "
<< toroidalResonance << "," << poloidalResonance << " ";
else if( type == FieldlineProperties::ISLANDS_WITHIN_ISLANDS )
else if( type == FieldlineProperties::ISLAND_WITH_SECONDARY_ISLANDS )
std::cerr << islands << " islands around "
<< islandGroups << " islandGroups with resonances: "
<< toroidalResonance << "," << poloidalResonance << " ";
else if( type == FieldlineProperties::ISLAND_AMBIGUOUS_AXIS )
{
std::cerr << islands << " island chain with an ambiguous axis: "
// << toroidalResonance << "," << poloidalResonance << " ";
<< toroidalWinding << "," << poloidalWindingP << " ("
<< (float) toroidalWinding / (float) poloidalWindingP << ") ";
}
else if( type == FieldlineProperties::CHAOTIC )
std::cerr << "chaotic ";
......@@ -819,12 +825,28 @@ avtPoincareFilter::CreatePoincareOutput( avtDataTree *dt,
<< std::endl;
if( (type == FieldlineProperties::ISLAND_CHAIN ||
type == FieldlineProperties::ISLANDS_WITHIN_ISLANDS) &&
type == FieldlineProperties::ISLAND_WITH_SECONDARY_ISLANDS ||
type == FieldlineProperties::ISLAND_AMBIGUOUS_AXIS) &&
toroidalWinding != poloidalWinding &&
islands != toroidalWinding )
std::cerr << "WARNING - The island count does not match the toroidalWinding count" << std::endl;
}
if( toroidalWinding == poloidalWinding )
{
if( type == FieldlineProperties::ISLAND_AMBIGUOUS_AXIS )
{
poloidalWinding = poloidalWindingP;
islands = 0;
}
else //if( type != FieldlineProperties::ISLAND_AMBIGUOUS_AXIS )
{
toroidalWinding = poloidalWinding = 1;
windingGroupOffset = 0;
}
}
// If toroidal winding is zero, skip it.
if( type == FieldlineProperties::CHAOTIC )
{
......@@ -1137,6 +1159,33 @@ avtPoincareFilter::CreatePoincareOutput( avtDataTree *dt,
}
}
}
else if( type == FieldlineProperties::ISLAND_AMBIGUOUS_AXIS )
{
if( overlaps != 0 )
{
if( showLines )
nnodes = 2;
bool tmpPoints = showPoints;
// Loop through each island.
for( unsigned int j=0; j<toroidalWinding; j++ )
{
// Erase all of the overlapping points.
puncturePts[p][j].erase( puncturePts[p][j].begin()+nnodes,
puncturePts[p][j].end() );
if( showLines )
{
unsigned int n = (j+windingGroupOffset) % toroidalWinding;
puncturePts[p][j][1] = puncturePts[p][j][0] +
0.9 * (puncturePts[p][n][0] - puncturePts[p][j][0]);
}
}
}
}
bool VALID = true;
......@@ -1171,7 +1220,7 @@ avtPoincareFilter::CreatePoincareOutput( avtDataTree *dt,
if( !showIslands ||
(showIslands &&
(type == FieldlineProperties::ISLAND_CHAIN ||
type == FieldlineProperties::ISLANDS_WITHIN_ISLANDS)) )
type == FieldlineProperties::ISLAND_WITH_SECONDARY_ISLANDS)) )
{
double color_value;
......@@ -1199,14 +1248,14 @@ avtPoincareFilter::CreatePoincareOutput( avtDataTree *dt,
}
else if( dataValue == DATA_SafetyFactorQ_NotP )
{
if( poloidalWinding == poloidalWindingP )
if( type != FieldlineProperties::ISLAND_AMBIGUOUS_AXIS )
color_value = (double) toroidalWinding / (double) poloidalWinding;
else
continue;
}
else if( dataValue == DATA_SafetyFactorP_NotQ )
{
if( poloidalWindingP != poloidalWinding )
if( type == FieldlineProperties::ISLAND_AMBIGUOUS_AXIS )
color_value = (double) toroidalWinding / (double) poloidalWindingP;
else
continue;
......@@ -1244,7 +1293,7 @@ avtPoincareFilter::CreatePoincareOutput( avtDataTree *dt,
windingGroupOffset,
dataValue, color_value );
}
else if( 0 && type == FieldlineProperties::ISLANDS_WITHIN_ISLANDS )
else if( 0 && type == FieldlineProperties::ISLAND_WITH_SECONDARY_ISLANDS )
{
drawIrrationalCurve( dt, puncturePts, nnodes, islands,
windingGroupOffset,
......@@ -1263,7 +1312,8 @@ avtPoincareFilter::CreatePoincareOutput( avtDataTree *dt,
if( showOPoints &&
(type == FieldlineProperties::ISLAND_CHAIN ||
type == FieldlineProperties::ISLANDS_WITHIN_ISLANDS) )
type == FieldlineProperties::ISLAND_WITH_SECONDARY_ISLANDS ||
type == FieldlineProperties::ISLAND_AMBIGUOUS_AXIS) )
{
drawPoints( dt, OPoints );
}
......@@ -1886,7 +1936,7 @@ avtPoincareFilter::drawSurface( avtDataTree *dt,
scalars->Allocate(dims[0]*dims[1]);
float *points_ptr = (float *) points->GetVoidPointer(0);
// Determine if the winding group order matches the point
// ordering. This is only needed when building surfaces.
Vector intra = nodes[0][ 0][1] - nodes[0][0][0];
......@@ -1919,7 +1969,7 @@ avtPoincareFilter::drawSurface( avtDataTree *dt,
{
k = j;
}
unsigned int jj = nplanes * j + p;
if( color == DATA_PlaneOrder )
......
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