avtPoincareFilter.C 161 KB
Newer Older
pugmire's avatar
pugmire committed
1 2
/*****************************************************************************
*
brugger's avatar
 
brugger committed
3
* Copyright (c) 2000 - 2013, Lawrence Livermore National Security, LLC
pugmire's avatar
pugmire committed
4
* Produced at the Lawrence Livermore National Laboratory
brugger's avatar
 
brugger committed
5
* LLNL-CODE-442911
pugmire's avatar
pugmire committed
6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
* All rights reserved.
*
* This file is  part of VisIt. For  details, see https://visit.llnl.gov/.  The
* full copyright notice is contained in the file COPYRIGHT located at the root
* of the VisIt distribution or at http://www.llnl.gov/visit/copyright.html.
*
* Redistribution  and  use  in  source  and  binary  forms,  with  or  without
* modification, are permitted provided that the following conditions are met:
*
*  - Redistributions of  source code must  retain the above  copyright notice,
*    this list of conditions and the disclaimer below.
*  - Redistributions in binary form must reproduce the above copyright notice,
*    this  list of  conditions  and  the  disclaimer (as noted below)  in  the
*    documentation and/or other materials provided with the distribution.
*  - Neither the name of  the LLNS/LLNL nor the names of  its contributors may
*    be used to endorse or promote products derived from this software without
*    specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT  HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR  PURPOSE
* ARE  DISCLAIMED. IN  NO EVENT  SHALL LAWRENCE  LIVERMORE NATIONAL  SECURITY,
* LLC, THE  U.S.  DEPARTMENT OF  ENERGY  OR  CONTRIBUTORS BE  LIABLE  FOR  ANY
* DIRECT,  INDIRECT,   INCIDENTAL,   SPECIAL,   EXEMPLARY,  OR   CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT  LIMITED TO, PROCUREMENT OF  SUBSTITUTE GOODS OR
* SERVICES; LOSS OF  USE, DATA, OR PROFITS; OR  BUSINESS INTERRUPTION) HOWEVER
* CAUSED  AND  ON  ANY  THEORY  OF  LIABILITY,  WHETHER  IN  CONTRACT,  STRICT
* LIABILITY, OR TORT  (INCLUDING NEGLIGENCE OR OTHERWISE)  ARISING IN ANY  WAY
* OUT OF THE  USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
* DAMAGE.
*
*****************************************************************************/

// ************************************************************************* //
allens's avatar
allens committed
40
//                             avtPoincareFilter.C                           //
pugmire's avatar
pugmire committed
41 42
// ************************************************************************* //

43 44 45

#define RATIONAL_SURFACE

pugmire's avatar
pugmire committed
46 47
#include <avtPoincareFilter.h>

allens's avatar
allens committed
48 49
#include <avtCallback.h>

50 51 52 53
#include <vtkAppendPolyData.h>
#include <vtkCellArray.h>
#include <vtkCellData.h>
#include <vtkCleanPolyData.h>
pugmire's avatar
pugmire committed
54
#include <vtkDataSet.h>
55 56
#include <vtkFloatArray.h>
#include <vtkPointData.h>
pugmire's avatar
pugmire committed
57
#include <vtkPolyData.h>
58
#include <vtkPolyLine.h>
59
#include <vtkQuad.h>
60 61 62
#include <vtkSlicer.h>
#include <vtkSphereSource.h>
#include <vtkTubeFilter.h>
63
#include <vtkUnstructuredGrid.h>
pugmire's avatar
pugmire committed
64

65 66
#include <avtDatasetExaminer.h>
#include <avtExtents.h>
67
#include <avtPoincareIC.h>
68 69
#include <utility>

allens's avatar
allens committed
70 71
#include <sys/stat.h>

allens's avatar
allens committed
72
#include "FieldlineAnalyzerLib.h"
73

allens's avatar
allens committed
74 75 76 77
#ifdef RATIONAL_SURFACE
#include "RationalSurfaceLib.h"
#endif

allens's avatar
allens committed
78
#ifdef STRAIGHTLINE_SKELETON
79
#include "skelet.h"
allens's avatar
allens committed
80
#endif
81

82
#define SIGN(x) ((x) < 0.0 ? -1 : 1)
pugmire's avatar
pugmire committed
83

allens's avatar
allens committed
84
static const int DATA_None = 0;
allens's avatar
allens committed
85 86 87 88 89 90 91 92 93 94 95 96 97
static const int DATA_SafetyFactorQ = 1;
static const int DATA_SafetyFactorP = 2;
static const int DATA_SafetyFactorQ_NotP = 3;
static const int DATA_SafetyFactorP_NotQ = 4;
static const int DATA_ToroidalWindings = 5;
static const int DATA_PoloidalWindingsQ = 6;
static const int DATA_PoloidalWindingsP = 7;
static const int DATA_FieldlineOrder = 8;
static const int DATA_PointOrder = 9;
static const int DATA_PlaneOrder = 10;
static const int DATA_WindingGroupOrder = 11;
static const int DATA_WindingPointOrder = 12;
static const int DATA_WindingPointOrderModulo = 13;
pugmire's avatar
pugmire committed
98

allens's avatar
allens committed
99 100 101 102 103 104 105 106
// ****************************************************************************
//  Method: CreateSphere
//
//  Programmer:
//  Creation:   Tue Oct 7 09:02:52 PDT 2008
//
//  Modifications:
//
pugmire's avatar
pugmire committed
107 108
//   Dave Pugmire, Thu Jul  1 13:55:28 EDT 2010
//   Create a vertex instead of a sphere.
109 110 111 112
//   
//    Dave Pugmire, Mon Jul 12 15:34:29 EDT 2010
//    Remove rad argument.
//    
allens's avatar
allens committed
113 114 115
// ****************************************************************************

static vtkPolyData *
116
CreateSphere(float val, double p[3])
allens's avatar
allens committed
117
{
pugmire's avatar
pugmire committed
118 119 120 121 122 123 124 125 126 127 128 129
    vtkPoints *pt = vtkPoints::New();
    pt->SetNumberOfPoints(1);
    pt->SetPoint(0, p[0], p[1], p[2]);
    
    vtkPolyData *point = vtkPolyData::New();
    point->SetPoints(pt);
    pt->Delete();

    vtkIdType ids[1] = {0};
    point->Allocate(1);
    point->InsertNextCell(VTK_VERTEX, 1, ids);

allens's avatar
allens committed
130 131
    vtkFloatArray *arr = vtkFloatArray::New();
    arr->SetName("colorVar");
pugmire's avatar
pugmire committed
132 133 134
    arr->SetNumberOfTuples(1);
    arr->SetTuple1(0, val);
    point->GetPointData()->SetScalars(arr);
allens's avatar
allens committed
135 136
    arr->Delete();

pugmire's avatar
pugmire committed
137
    return point;
allens's avatar
allens committed
138 139 140
}


pugmire's avatar
pugmire committed
141 142 143
// ****************************************************************************
//  Method: avtPoincareFilter constructor
//
144
//  Programmer: Dave Pugmire -- generated by xml2avt
pugmire's avatar
pugmire committed
145 146
//  Creation:   Tue Oct 7 09:02:52 PDT 2008
//
147 148 149
//  Modifications:
//
//    Dave Pugmire, Wed Feb 25 09:52:11 EST 2009
150 151
//    Add terminate by steps, add AdamsBashforth solver,
//    Allen Sanderson's new code.
152
//
pugmire's avatar
pugmire committed
153
//    Dave Pugmire, Fri Apr 17 11:32:40 EDT 2009
154
//    Add variables for dataValue var.
pugmire's avatar
pugmire committed
155
//
pugmire's avatar
pugmire committed
156
//    Dave Pugmire, Tue Apr 28 09:26:06 EDT 2009
157
//    Changed color to dataValue
pugmire's avatar
pugmire committed
158
//
pugmire's avatar
pugmire committed
159
//    Dave Pugmire, Wed May 27 15:03:42 EDT 2009
160
//    Initialize points.
pugmire's avatar
pugmire committed
161
//
pugmire's avatar
pugmire committed
162
//    Dave Pugmire, Tue Aug 18 09:10:49 EDT 2009
163
//    Add ability to restart fieldline integration.
pugmire's avatar
pugmire committed
164
//
pugmire's avatar
pugmire committed
165 166
// ****************************************************************************

allens's avatar
allens committed
167
avtPoincareFilter::avtPoincareFilter() :
168 169
    maximumToroidalWinding( 0 ),
    overrideToroidalWinding( 0 ),
allens's avatar
allens committed
170
    overridePoloidalWinding( 0 ),
171
    windingPairConfidence( 0.90 ),
allens's avatar
allens committed
172
    rationalSurfaceFactor( 0.10 ),
173
    adjust_plane(-1),
174
    overlaps(1),
175

176
    is_curvemesh(1),
177
    dataValue(DATA_SafetyFactorQ),
178

allens's avatar
allens committed
179 180
    showRationalSurfaces( false ),
    rationalSurfaceMaxIterations(2),
181
    showOPoints( false ),
allens's avatar
allens committed
182 183
    OPointMaxIterations(2),
    XPointMaxIterations(2),
allens's avatar
allens committed
184
    performOLineAnalysis( false ),
allens's avatar
allens committed
185
    OLineToroidalWinding( 1 ),
allens's avatar
allens committed
186
    OLineAxisFileName(""),
187 188 189
    showIslands( false ),
    showLines( true ),
    showPoints( false ),
allens's avatar
allens committed
190 191
    summaryFlag( true ),
    verboseFlag( false ),
192
    pointScale(1)
pugmire's avatar
pugmire committed
193
{
194 195
    planes.resize(1);
    planes[0] = 0;
196 197
    intersectObj = NULL;
    maxIntersections = 0;
pugmire's avatar
pugmire committed
198 199 200 201 202 203
}


// ****************************************************************************
//  Method: avtPoincareFilter destructor
//
204
//  Programmer: Dave Pugmire -- generated by xml2avt
pugmire's avatar
pugmire committed
205 206
//  Creation:   Tue Oct 7 09:02:52 PDT 2008
//
pugmire's avatar
pugmire committed
207 208
//  Modifications:
//
pugmire's avatar
pugmire committed
209 210 211 212
// ****************************************************************************

avtPoincareFilter::~avtPoincareFilter()
{
213 214
    if (intersectObj)
        intersectObj->Delete();
pugmire's avatar
pugmire committed
215 216
}

217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248
// ****************************************************************************
//  Method: avtPoincareFilter::PreExecute
//
//  Purpose:
//      Get the current spatial extents if necessary.
//
//  Programmer: Dave Pugmire
//  Creation:   Fri Nov  7 13:01:47 EST 2008
//
//  Modifications:
//
//
// ****************************************************************************

void
avtPoincareFilter::PreExecute(void)
{
    avtStreamlineFilter::PreExecute();
}


// ****************************************************************************
//  Method: avtPoincareFilter::PostExecute
//
//  Purpose:
//      Get the current spatial extents if necessary.
//
//  Programmer: Dave Pugmire
//  Creation:   Fri Nov  7 13:01:47 EST 2008
//
//  Modifications:
//
249 250
//    Hank Childs, Thu Aug 26 13:47:30 PDT 2010
//    Change extents names.
251 252
//
// ****************************************************************************
pugmire's avatar
pugmire committed
253

254 255 256 257
void
avtPoincareFilter::PostExecute(void)
{
    avtStreamlineFilter::PostExecute();
258
    
259 260 261
    double range[2];
    avtDataset_p ds = GetTypedOutput();
    avtDatasetExaminer::GetDataExtents(ds, range, "colorVar");
pugmire's avatar
pugmire committed
262

263
    avtExtents *e;
264
    e = GetOutput()->GetInfo().GetAttributes().GetThisProcsOriginalDataExtents();
265
    e->Merge(range);
266
    e = GetOutput()->GetInfo().GetAttributes().GetThisProcsActualDataExtents();
267 268
    e->Merge(range);
}
pugmire's avatar
pugmire committed
269

270 271 272 273 274 275 276 277
// ****************************************************************************
//  Method: avtPoincareFilter::CreateIntegralCurve
//
//  Purpose:
//      Create an integral curve and set its properties.
//
//  Programmer: Christoph Garth
//  Creation:   Thu July 15, 2010
278 279 280 281 282 283
//
//  Modifications:
//
//    Hank Childs, Fri Oct  8 23:30:27 PDT 2010
//    Create PoincareICs, not StateRecorderICs.
//
284 285 286 287 288 289
// ****************************************************************************

avtIntegralCurve *
avtPoincareFilter::CreateIntegralCurve( const avtIVPSolver* model,
                                        avtIntegralCurve::Direction dir,
                                        const double& t_start,
290 291 292
                                        const avtVector &p_start,
                                        const avtVector &v_start,
                                        long ID ) 
293 294 295 296
{
    // need at least these three attributes
    unsigned char attr = avtStateRecorderIntegralCurve::SAMPLE_POSITION;

297
    avtPoincareIC *rv = new avtPoincareIC( attr, model, dir, 
298
                                           t_start, p_start, v_start, ID );
299 300

    if (intersectObj)
301
        rv->SetIntersectionCriteria(intersectObj, maxIntersections);
302 303 304 305

    return rv;
}

306
// ****************************************************************************
307
//  Method: avtPoincareFilter::GetFieldlinePoints
308 309
//
//  Purpose:
310
//      Gets the points from the fieldline and changes them in to a Vector.
311 312 313 314 315 316 317 318 319
//
//  Programmer: Dave Pugmire
//  Creation:   Tue Dec  23 12:51:29 EST 2008
//
//  Modifications:
//
// ****************************************************************************

void
allens's avatar
allens committed
320
avtPoincareFilter::GetIntegralCurvePoints(std::vector<avtIntegralCurve *> &ics)
321
{
322
    for ( int i=0; i<ics.size(); ++i )
pugmire's avatar
pugmire committed
323
    {
324
        avtPoincareIC * poincare_ic = (avtPoincareIC *) ics[i];
allens's avatar
allens committed
325

326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349
        // Only move the points needed over to the storage.
        if( poincare_ic->points.size() < poincare_ic->GetNumberOfSamples() )
        {
          unsigned int start = poincare_ic->points.size();
          unsigned int stop  = poincare_ic->GetNumberOfSamples();

          // Get all of the points from the fieldline which are stored
          // as an array and move them into a vector for easier
          // manipulation by the analysis code.
          poincare_ic->points.resize( poincare_ic->GetNumberOfSamples() );

          for( size_t p=start; p<stop; ++p )
            poincare_ic->points[p] = poincare_ic->GetSample( p ).position;
        }

        // If the analysis asked for more points but did not get any
        // then the integration failed. As such, terminate the
        // integration and analysis.
        else if( poincare_ic->properties.analysisState ==
                 FieldlineProperties::ADDING_POINTS &&
                 
                 poincare_ic->points.size() ==
                 poincare_ic->GetNumberOfSamples() )
        {
350 351
            poincare_ic->status.SetTerminationMet();
            poincare_ic->properties.analysisState =
352
            FieldlineProperties::TERMINATED;
allens's avatar
allens committed
353

354 355 356 357
//           std::cerr << "Terminated integration for Fieldline: id = "
//                     << poincare_ic->id << "  "
//                  << std::endl;
        }
allens's avatar
allens committed
358
    }
359 360 361 362 363 364 365 366 367 368 369
}

// ****************************************************************************
//  Method: avtPoincareFilter::Execute
//
//  Purpose:
//      Calculate poincare points.
//
//  Programmer: Dave Pugmire -- generated by xml2avt
//  Creation:   Tue Oct 7 09:02:52 PDT 2008
//
370 371 372 373
//  Modifications:
//    Dave Pugmire (for Allen Sanderson), Wed Feb 25 09:52:11 EST 2009
//    Add terminate by steps, add AdamsBashforth solver, Allen Sanderson's new code.
//
pugmire's avatar
pugmire committed
374
//    Dave Pugmire, Wed May 27 15:03:42 EDT 2009
375
//    Re-organization. GetFieldlinePoints removed.
pugmire's avatar
pugmire committed
376
//
pugmire's avatar
pugmire committed
377
//    Dave Pugmire, Tue Aug 18 09:10:49 EDT 2009
378
//    Add ability to restart fieldline integration.
pugmire's avatar
pugmire committed
379
//
380 381 382 383 384
// ****************************************************************************

void
avtPoincareFilter::Execute()
{
allens's avatar
allens committed
385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426
    if( performOLineAnalysis )
    {
        struct stat fileAtt;

        //Use the stat function to get the file information
        if (stat(OLineAxisFileName.c_str(), &fileAtt) != 0)
        {
          std::string msg("Trying to perform O-line analysis but the O-line axis file is not valid.");

          avtCallback::IssueWarning(msg.c_str());
          EXCEPTION1(ImproperUseException, msg);
        }
        else
        {
          FILE *fp = fopen( OLineAxisFileName.c_str(), "r" );

          if( fp == NULL )
          {
            std::string msg("Trying to perform O-line analysis but the O-line axis file can not be openned.");
            
            avtCallback::IssueWarning(msg.c_str());
            EXCEPTION1(ImproperUseException, msg);
          }

          while( !feof(fp) )
          {
            Point nextPt;

            if( fscanf( fp, "%lf %lf %lf", &nextPt.x, &nextPt.y, &nextPt.z ) != 3 &&
                feof(fp) == 0)
            {
              std::string msg("Trying to perform O-line analysis but the O-line axis file can not be read.");
              
              avtCallback::IssueWarning(msg.c_str());
              EXCEPTION1(ImproperUseException, msg);
            }
          }

          fclose (fp);
        }
    }

427
    avtStreamlineFilter::Execute();
428

allens's avatar
allens committed
429
    std::vector<avtIntegralCurve *> ics;
430 431
    GetTerminatedIntegralCurves(ics);

allens's avatar
allens committed
432 433 434
    avtDataTree *dt = new avtDataTree();
    
    CreatePoincareOutput( dt, ics );
allens's avatar
allens committed
435 436 437
#ifdef RATIONAL_SURFACE
    CreateRationalOutput( dt, ics );
#endif
pugmire's avatar
pugmire committed
438
    SetOutputDataTree(dt);
pugmire's avatar
pugmire committed
439 440
}

allens's avatar
allens committed
441

pugmire's avatar
pugmire committed
442 443 444 445 446 447 448 449 450 451 452
// ****************************************************************************
//  Method: avtPoincareFilter::ContinueExecute
//
//  Purpose:
//      See if execution needs to continue.
//
//  Programmer: Dave Pugmire
//  Creation:   Mon Aug 17 08:30:06 EDT 2009
//
//  Modifications:
//
453
//    Hank Childs, Sun Jun  6 11:53:33 CDT 2010
454
//    Use new names that have integral curves instead of fieldlines.
455
//
pugmire's avatar
pugmire committed
456 457 458 459 460
// ****************************************************************************

bool
avtPoincareFilter::ContinueExecute()
{
allens's avatar
allens committed
461
    debug5 << "Continue execute " << std::endl;
462

allens's avatar
allens committed
463
    std::vector<avtIntegralCurve *> ics;
pugmire's avatar
pugmire committed
464
    
465 466
    GetTerminatedIntegralCurves(ics);
    GetIntegralCurvePoints(ics);
pugmire's avatar
pugmire committed
467

468
    if (analysis && (! ClassifyFieldlines(ics) || ! ClassifyRationals(ics)))
pugmire's avatar
pugmire committed
469
    {
allens's avatar
allens committed
470
      std::vector< int > ids_to_delete;
allens's avatar
allens committed
471

allens's avatar
allens committed
472
      // Because points are added the size will change so get the
473
      // inital size so that new seeds are not processed.
allens's avatar
allens committed
474 475 476
      unsigned int nics = ics.size();

      for ( int i=0; i<nics; ++i )
allens's avatar
allens committed
477
      {
478
        avtPoincareIC * poincare_ic = (avtPoincareIC *) ics[i];
479
        FieldlineProperties &properties = poincare_ic->properties;
480

allens's avatar
allens committed
481
#ifdef STRAIGHTLINE_SKELETON
allens's avatar
allens committed
482
        // For Island Chains add in the O Points.
allens's avatar
allens committed
483
        if( showOPoints )
allens's avatar
allens committed
484
        {
allens's avatar
allens committed
485 486 487 488
          // Are O Points present?
          if( properties.type & FieldlineProperties::ISLAND_CHAIN &&
              properties.analysisState == FieldlineProperties::ADD_O_POINTS &&

489 490
              (properties.searchState == FieldlineProperties::ISLAND_O_POINT ||
               properties.searchState == FieldlineProperties::NO_SEARCH) )
allens's avatar
allens committed
491
          {
allens's avatar
allens committed
492 493 494 495 496 497 498 499 500 501
            // Change the state of the properties to complete now that
            // the seed point has been stripped off.
            poincare_ic->properties.analysisState =
              FieldlineProperties::COMPLETED;
            
            if(verboseFlag )
              std::cerr << "Have island seed  " << properties.seedPoints[0]
                        << std::endl;
            
            if( properties.iteration < OPointMaxIterations )
allens's avatar
allens committed
502
            {
allens's avatar
allens committed
503 504 505
              std::vector<avtIntegralCurve *> new_ics;
              avtVector vec(0,0,0);
              
allens's avatar
allens committed
506
              if(verboseFlag )
allens's avatar
allens committed
507 508
                std::cerr << "Adding island O point seed  "
                          << properties.seedPoints[0]
allens's avatar
allens committed
509
                          << std::endl;
allens's avatar
allens committed
510 511 512 513 514 515 516 517 518 519 520 521 522 523 524
              
              AddSeedPoint( properties.seedPoints[0], vec, new_ics );
              
              for( unsigned int j=0; j<new_ics.size(); ++j )
              {
                if(verboseFlag )
                  std::cerr << "New island O point seed ids "
                            << new_ics[j]->id << "  ";

                avtPoincareIC* seed_poincare_ic =
                  (avtPoincareIC *) new_ics[j];

                // Transfer and update properties.
                seed_poincare_ic->properties = properties;
              
525 526
                seed_poincare_ic->properties.searchState =
                  FieldlineProperties::ISLAND_O_POINT;
allens's avatar
allens committed
527
                seed_poincare_ic->properties.analysisState =
528 529 530
                  FieldlineProperties::UNKNOWN_ANALYSIS;
                seed_poincare_ic->properties.source =
                  properties.type;
allens's avatar
allens committed
531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546
              
                seed_poincare_ic->properties.iteration =
                  properties.iteration + 1;
              }

              if(verboseFlag )
                std::cerr << std::endl;
              
              // The source was an island_chain which meant the seed was
              // an intermediate seed so delete it.

              // Note only delete the seed if another seed replaces
              // it. If past the maximum iterations the seed will
              // not be deleted.
              if( properties.source & FieldlineProperties::ISLAND_CHAIN )
              {
allens's avatar
allens committed
547
                if( verboseFlag )
548
                  std::cerr << "O Point search deleting O Point seed "
allens's avatar
allens committed
549 550 551 552 553 554 555
                            << poincare_ic->id << std::endl;
                
                ids_to_delete.push_back( poincare_ic->id );
              }
            }
          }
              
556 557
          // Landed on an O-point directly. Can not search for the
          // boundary as there is no value for properties.searchDelta.
allens's avatar
allens committed
558
          else if( properties.type == FieldlineProperties::O_POINT &&
559 560
                   properties.analysisState == FieldlineProperties::ADD_BOUNDARY_POINT &&
                   properties.searchState == FieldlineProperties::NO_SEARCH )
allens's avatar
allens committed
561 562 563 564 565
          {
            // Change the state of the properties to complete now that
            // the seed point has been stripped off.
            poincare_ic->properties.analysisState =
              FieldlineProperties::COMPLETED;
566
          }
allens's avatar
allens committed
567

568 569 570
          // Is a boundary seed point present?
          else if( properties.type == FieldlineProperties::O_POINT &&
                   properties.analysisState == FieldlineProperties::ADD_BOUNDARY_POINT &&
allens's avatar
allens committed
571

572 573 574 575 576 577 578
                    (properties.searchState == FieldlineProperties::ISLAND_O_POINT ||
                     properties.searchState == FieldlineProperties::ISLAND_BOUNDARY_SEARCH) )
          {
            // Change the state of the properties to complete now that
            // the seed point has been stripped off.
            poincare_ic->properties.analysisState =
              FieldlineProperties::COMPLETED;
allens's avatar
allens committed
579 580

            properties.searchIncrement = 1.0;
581
            properties.pastFirstSearchFailure = false;
allens's avatar
allens committed
582

583
            if( properties.iteration < OPointMaxIterations )
allens's avatar
allens committed
584 585 586 587
            {
              // First time through the loop.
              if( properties.searchState == FieldlineProperties::ISLAND_O_POINT )
                properties.searchMagnitude = properties.searchIncrement;
588

allens's avatar
allens committed
589
              // Ended up back up at the O Point so try again.
590
              else if( properties.searchState == FieldlineProperties::ISLAND_BOUNDARY_SEARCH )
allens's avatar
allens committed
591 592 593 594 595 596 597 598 599
                properties.searchMagnitude += properties.searchIncrement;

              properties.baseToroidalWinding = properties.toroidalWinding;
              properties.basePoloidalWinding = properties.poloidalWinding;
              
              avtVector seed = properties.lastSeedPoint +
                properties.searchNormal *
                properties.searchMagnitude * properties.searchDelta;
              
600 601 602 603 604 605 606 607 608 609
              std::cerr << "LINE " << __LINE__ << "  "
                        << properties.iteration << "  "
                        << properties.pastFirstSearchFailure << "  "
                        << properties.searchIncrement << "  "
                        << properties.searchMagnitude << "  "
                        << seed << "  "
                        << properties.searchNormal << "  "
                        << properties.searchDelta << "  "
                        << std::endl;

allens's avatar
allens committed
610 611 612 613
              std::vector<avtIntegralCurve *> new_ics;
              avtVector vec(0,0,0);
              
              if(verboseFlag )
614
                std::cerr << "Have island boundary seed  " << seed << std::endl;
allens's avatar
allens committed
615 616 617 618 619 620
              
              AddSeedPoint( seed, vec, new_ics );
              
              for( unsigned int j=0; j<new_ics.size(); ++j )
              {
                if(verboseFlag )
621 622
                  std::cerr << "LINE " << __LINE__
                            << " New island boundary seed ids "
allens's avatar
allens committed
623 624 625 626 627 628 629 630
                            << new_ics[j]->id << "  ";
                
                avtPoincareIC* seed_poincare_ic =
                  (avtPoincareIC *) new_ics[j];
                
                // Transfer and update properties.
                seed_poincare_ic->properties = properties;
                
631 632
                seed_poincare_ic->properties.searchState =
                  FieldlineProperties::ISLAND_BOUNDARY_SEARCH;
allens's avatar
allens committed
633
                seed_poincare_ic->properties.analysisState =
634
                  FieldlineProperties::UNKNOWN_ANALYSIS;
allens's avatar
allens committed
635 636
                
                seed_poincare_ic->properties.source = properties.type;
637 638
                // Save the seed point curve.
                seed_poincare_ic->properties.parentOPointIC = poincare_ic;
allens's avatar
allens committed
639
              
640
                // First time through the loop.
allens's avatar
allens committed
641
                if( properties.searchState == FieldlineProperties::ISLAND_O_POINT )
642
                {
allens's avatar
allens committed
643
                  seed_poincare_ic->properties.iteration = 0;
644
                }
allens's avatar
allens committed
645

646 647 648
                // Ended up back up at the O Point so try again with a
                // new seed.
                else if( properties.searchState == FieldlineProperties::ISLAND_BOUNDARY_SEARCH )
allens's avatar
allens committed
649
                {
allens's avatar
allens committed
650 651 652
                  seed_poincare_ic->properties.iteration =
                    properties.iteration + 1;

allens's avatar
allens committed
653
                  if( verboseFlag )
654
                    std::cerr << "Island boundary search deleting O Point seed "
allens's avatar
allens committed
655 656 657 658 659
                              << poincare_ic->id << std::endl;
                
                  ids_to_delete.push_back( poincare_ic->id );
                }
              }
allens's avatar
allens committed
660 661 662 663 664 665

              if(verboseFlag )
                std::cerr << std::endl;
            }
          }

666
          // Boundary surfaces
allens's avatar
allens committed
667
          else if( ( (properties.type & FieldlineProperties::ISLAND_CHAIN &&
668 669 670 671 672 673 674 675
                      properties.analysisState == FieldlineProperties::ADD_BOUNDARY_POINT) || 

                     (properties.type & FieldlineProperties::FLUX_SURFACE &&
                      properties.analysisState == FieldlineProperties::COMPLETED) || 
                     
                     (properties.analysisState == FieldlineProperties::TERMINATED) ) &&

                   properties.searchState == FieldlineProperties::ISLAND_BOUNDARY_SEARCH )
allens's avatar
allens committed
676
          {
677 678
            bool terminated = properties.analysisState == FieldlineProperties::TERMINATED;

allens's avatar
allens committed
679 680 681 682
            // Change the state of the properties to complete.
            poincare_ic->properties.analysisState =
              FieldlineProperties::COMPLETED;

683 684
            if( properties.iteration <
                (properties.pastFirstSearchFailure ? 1 : 10) * OPointMaxIterations )
allens's avatar
allens committed
685
            {
686 687 688 689
              // If a flux surface is found or the analysis terminates
              // go back a step and then half the increment. 
              if( (properties.type & FieldlineProperties::FLUX_SURFACE) ||
                  terminated )
allens's avatar
allens committed
690 691 692 693 694 695 696 697 698
              {
                if( properties.pastFirstSearchFailure == false )
                {
                  properties.pastFirstSearchFailure = true;
                  properties.iteration = 0;
                }

                properties.searchMagnitude -= properties.searchIncrement;
              }
699 700 701

              if( properties.pastFirstSearchFailure )
                properties.searchIncrement /= 2.0;
allens's avatar
allens committed
702 703 704
              
              // If about to end do not increment so to be assured
              // that an island is found.
705 706
              if( properties.iteration + 1 <
                  (properties.pastFirstSearchFailure ? 1 : 10) * OPointMaxIterations )
allens's avatar
allens committed
707 708 709 710 711 712
                properties.searchMagnitude += properties.searchIncrement;

              avtVector seed = properties.lastSeedPoint +
                properties.searchMagnitude * properties.searchDelta *
                properties.searchNormal;
            
713
              std::cerr << "LINE " << __LINE__ << "  "
allens's avatar
allens committed
714
                        << properties.iteration << "  "
715 716 717
                        << properties.pastFirstSearchFailure << "  "
                        << properties.searchIncrement << "  "
                        << properties.searchMagnitude << "  "
allens's avatar
allens committed
718 719 720 721 722 723 724 725 726
                        << seed << "  "
                        << properties.searchNormal << "  "
                        << properties.searchDelta << "  "
                        << std::endl;

              std::vector<avtIntegralCurve *> new_ics;
              avtVector vec(0,0,0);
              
              if(verboseFlag )
727 728
                std::cerr << "LINE " << __LINE__
                          << " Have additional island boundary seed  " << seed << std::endl;
allens's avatar
allens committed
729 730 731 732 733 734
                          
              AddSeedPoint( seed, vec, new_ics );
              
              for( unsigned int j=0; j<new_ics.size(); ++j )
              {
                if(verboseFlag )
735
                  std::cerr << "New island boundary seed ids "
allens's avatar
allens committed
736 737 738 739 740 741 742 743 744
                            << new_ics[j]->id << "  ";
              
                avtPoincareIC* seed_poincare_ic =
                  (avtPoincareIC *) new_ics[j];
              
                // Transfer and update properties.
                seed_poincare_ic->properties = properties;
              
                seed_poincare_ic->properties.analysisState =
745
                  FieldlineProperties::UNKNOWN_ANALYSIS;
allens's avatar
allens committed
746 747 748 749 750 751 752
              
                seed_poincare_ic->properties.source = properties.type;
              
                seed_poincare_ic->properties.iteration =
                  properties.iteration + 1;
              
                seed_poincare_ic->properties.searchState =
753
                  FieldlineProperties::ISLAND_BOUNDARY_SEARCH;
allens's avatar
allens committed
754 755 756 757 758 759 760 761 762 763
              }

              if(verboseFlag )
                std::cerr << std::endl;

              // Note only delete the seed if another seed replaces
              // it. If past the maximum iterations the seed will
              // not be deleted.
//            if( properties.source & FieldlineProperties::ISLAND_CHAIN )
              {
allens's avatar
allens committed
764
                if( verboseFlag )
765
                  std::cerr << "Island boundary search deleting boundary seed "
allens's avatar
allens committed
766 767 768 769 770 771 772 773
                            << poincare_ic->id << std::endl;
                
                ids_to_delete.push_back( poincare_ic->id );
              }
            }
            else
            {
              poincare_ic->properties.searchState =
774
                FieldlineProperties::NO_SEARCH;
allens's avatar
allens committed
775

776
              std::cerr << "LINE " << __LINE__ << "  "
allens's avatar
allens committed
777 778 779 780 781 782 783 784 785
                        << properties.baseToroidalWinding << "  "
                        << properties.basePoloidalWinding << "  "
                        << properties.toroidalWinding << "  "
                        << properties.poloidalWinding << "  "
                        << std::endl;

              poincare_ic->properties.parentOPointIC->properties.childOPointIC = 
                poincare_ic;
            }
allens's avatar
allens committed
786 787
          }

allens's avatar
allens committed
788 789 790 791
          // Not an O or X point check to see if the source was from
          // an island chain. If so delete.
          else if( properties.type != FieldlineProperties::O_POINT &&
                   properties.type != FieldlineProperties::X_POINT )
allens's avatar
allens committed
792
          {
allens's avatar
allens committed
793 794 795
            // The source was an island_chain which meant the seed was
            // an intermediate seed that did make it into an O Point
            // so delete it.
allens's avatar
allens committed
796
            if( properties.source & FieldlineProperties::ISLAND_CHAIN &&
797
                properties.analysisState == FieldlineProperties::COMPLETED &&
allens's avatar
allens committed
798
                properties.searchState == FieldlineProperties::ISLAND_O_POINT )
allens's avatar
allens committed
799
            {
allens's avatar
allens committed
800
              if( verboseFlag )
801
                std::cerr << "Deleting an O Point seed that resulted in a surface "
allens's avatar
allens committed
802 803 804 805
                          << poincare_ic->id << std::endl;
              
              ids_to_delete.push_back( poincare_ic->id );
            }
allens's avatar
allens committed
806 807 808
          }
        }
#endif
allens's avatar
allens committed
809 810

#ifdef RATIONAL_SURFACE
allens's avatar
allens committed
811 812 813 814

        if( showRationalSurfaces )
        {

allens's avatar
allens committed
815 816 817
        /////////////////////////
        // Begin Rational Search
        /////////////////////////
818 819 820 821 822 823 824
        if ( properties.analysisMethod == FieldlineProperties::RATIONAL_SEARCH &&
             properties.type        == FieldlineProperties::RATIONAL &&
             properties.searchState == FieldlineProperties::ORIGINAL_RATIONAL )
          {       
            // Intentionally empty
          }
        else if( properties.analysisMethod  == FieldlineProperties::DEFAULT_METHOD &&
825 826 827
            properties.type            == FieldlineProperties::RATIONAL &&
            (properties.analysisState  == FieldlineProperties::COMPLETED ||
             properties.analysisState  == FieldlineProperties::TERMINATED) )
allens's avatar
allens committed
828
        {
829

830 831 832
          // Here we've caught a rational coming out of analysis
          //Set rational to original_rational, and rational_search and setup children vector
          SetupRational(poincare_ic);
833

834 835
          avtVector point1, point2;
          std::vector<avtVector> seeds = GetSeeds(poincare_ic, point1, point2, MAX_SEED_SPACING);
836

837 838 839 840 841 842 843 844 845 846 847 848 849 850
          for( unsigned int s=0; s<seeds.size(); ++s )
            {   
              if (2 <= RATIONAL_DEBUG)std::cerr << "LINE " << __LINE__ << " New Seed: " 
                                                << VectorToString(seeds[s]);
             
              std::vector<avtIntegralCurve *> new_ics;
              avtVector vec(0,0,0);
              
              seeds[s][1] = Z_OFFSET;
              AddSeedPoint( seeds[s], vec, new_ics );
              
              for( unsigned int j=0; j<new_ics.size(); ++j )
                {
                  avtPoincareIC *seed = (avtPoincareIC *) new_ics[j];
851

852 853 854 855 856 857
                  SetupNewSeed(seed, poincare_ic, seeds[s], point1, point2);              
                  
                  if (2 <= RATIONAL_DEBUG)std::cerr << " with ID: " << seed->id <<"\n";
                }
              if (2 <= RATIONAL_DEBUG)std::cerr << std::endl;
            }
allens's avatar
allens committed
858 859
        }

860 861 862 863 864 865 866
        //////////////////
        // Grab each seed 
        //////////////////
        else if( properties.analysisMethod == FieldlineProperties::RATIONAL_SEARCH &&
                 properties.searchState    == FieldlineProperties::SEARCHING_SEED &&
                 (properties.analysisState  == FieldlineProperties::COMPLETED ||
                  properties.analysisState  == FieldlineProperties::TERMINATED))
allens's avatar
allens committed
867
          {
868 869
          if (1 <= RATIONAL_DEBUG)std::cerr << "LINE: " << __LINE__ << "  "
                                            << "Found seed ID: " << poincare_ic->id <<",pt: "<< VectorToString(poincare_ic->points[0])<< std::endl;
870

871 872 873 874 875
          if(NO_MINIMIZATION)
            {
              poincare_ic->properties.searchState = FieldlineProperties::WAITING_SEED;
              
            }
allens's avatar
allens committed
876 877
          else
            {
878 879
              bool needToMinimize = NeedToMinimize(poincare_ic);
              if( !needToMinimize )
allens's avatar
allens committed
880
                {
881 882 883 884
                  if (1 <= RATIONAL_DEBUG)std::cerr<< "LINE: " << __LINE__
                                                   << " Got lucky, this seed ID: " <<poincare_ic->id
                                                   <<"  is already very good: " << FindMinimizationDistance(poincare_ic)<<std::endl;
                  properties.searchState = FieldlineProperties::WAITING_SEED;
allens's avatar
allens committed
885
                }
886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925
              else
                {          // Otherwise, we need to bracket the min before we can minimize to it
                  // May as well guess our seed is nearest the minimum, so make it b
                  avtVector zeroVec = avtVector(0,0,0);
              
                  avtVector newA, newB, newC; //newA isn't used, poincare_ic is A
                  if (false == PrepareToBracket(poincare_ic, newA, newB, newC))
                    {
                      poincare_ic->properties.searchState = FieldlineProperties::WAITING_SEED;
                      /* std::vector<avtPoincareIC *> *children = poincare_ic->src_rational_ic->properties.children;
                      if (2 <= RATIONAL_DEBUG)
                        std::cerr << "LINE " << __LINE__ << "  " << "Failed to prepare bracketing, Deleting ID: "
                                  <<poincare_ic->id << std::endl;
                      ids_to_delete.push_back(poincare_ic->id);
                      if(std::find(children->begin(), children->end(), poincare_ic) != children->end())
                      children->erase(std::remove(children->begin(), children->end(), poincare_ic));*/
                    }
                  else
                    {
                      std::vector<avtIntegralCurve *> new_ics;
                      newA[1] =Z_OFFSET;                      
                      AddSeedPoint( newB, zeroVec, new_ics );
                      avtPoincareIC *seed_b;          
                      for( unsigned int k=0; k<new_ics.size(); k++ )
                        {
                          seed_b = (avtPoincareIC *) new_ics[k];
                          SetupNewBracketB(seed_b, poincare_ic, newB);
                        }
                  
                      // Setup 'c' (use the new point just calculated)
                      std::vector<avtIntegralCurve *> new_ics_2;
                      newC[1]=Z_OFFSET;
                      AddSeedPoint( newC, zeroVec, new_ics_2 );
                      avtPoincareIC *seed_c;
                  
                      for( unsigned int k=0; k<new_ics_2.size(); k++ )
                        {
                          seed_c = (avtPoincareIC *) new_ics_2[k];
                          SetupNewBracketC(seed_c, poincare_ic, newC);
                        }
allens's avatar
allens committed
926

927 928 929
                      // We catch the b curve when it comes back to bracket the minimum
                      seed_b->a_IC = poincare_ic;
                      seed_b->c_IC = seed_c;
930

931
                    }     
allens's avatar
allens committed
932 933
                }
            }
934
        }      
935 936
        else if (properties.analysisMethod == FieldlineProperties::RATIONAL_SEARCH &&
                 properties.searchState == FieldlineProperties::MINIMIZING_A)
937 938 939 940 941
          {
            // Intentionally empty
            if (5 <= RATIONAL_DEBUG)
              cerr << "Minimizing A found" << std::endl;
          }
942 943
        else if (properties.analysisMethod == FieldlineProperties::RATIONAL_SEARCH &&
                 properties.searchState == FieldlineProperties::MINIMIZING_C)
944 945 946 947 948
          {
            // Intentionally empty
            if (5 <= RATIONAL_DEBUG)
              cerr << "Minimizing C found" << std::endl;
          }
949 950
        else if ( properties.analysisMethod == FieldlineProperties::RATIONAL_SEARCH &&
                  properties.searchState == FieldlineProperties::MINIMIZING_B &&
951 952
                  poincare_ic->a_IC != NULL && poincare_ic->c_IC != NULL &&
                  std::find(ids_to_delete.begin(), ids_to_delete.end(), poincare_ic->id) == ids_to_delete.end())
allens's avatar
allens committed
953 954 955
        {
          avtPoincareIC *_a = poincare_ic->a_IC;
          avtPoincareIC *_b = poincare_ic;
956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971
          avtPoincareIC *_c = poincare_ic->c_IC;
          avtPoincareIC *seed = poincare_ic->src_seed_ic;
          
          Vector xzplane(0,1,0);
          FieldlineLib fieldlib;
          std::vector<avtVector> xa_puncturePoints;
          std::vector<avtVector> xb_puncturePoints;
          std::vector<avtVector> xc_puncturePoints;       
          fieldlib.getPunctures(_a->points,xzplane,xa_puncturePoints);
          fieldlib.getPunctures(_b->points,xzplane,xb_puncturePoints);
          fieldlib.getPunctures(_c->points,xzplane,xc_puncturePoints);
          int xa_i = FindMinimizationIndex(_a);
          int xb_i = FindMinimizationIndex(_b);
          int xc_i = FindMinimizationIndex(_c);

          if (1 <= RATIONAL_DEBUG)
allens's avatar
allens committed
972
            {
973 974 975 976
              std::cerr <<"Line: "<<__LINE__<< " Found Bracketing\n"
                        <<"Line: "<<__LINE__<< " A ID: " <<_a->id <<", dist: "<< FindMinimizationDistance(_a) <<"\tpt: " <<VectorToString(xa_puncturePoints[xa_i])<< std::endl;
              std::cerr <<"Line: "<<__LINE__<< " B ID: " <<_b->id <<", dist: "<< FindMinimizationDistance(_b) <<"\tpt: " <<VectorToString(xb_puncturePoints[xb_i])<< std::endl;
              std::cerr <<"Line: "<<__LINE__<< " C ID: "        <<_c->id <<", dist: "<< FindMinimizationDistance(_c) <<"\tpt: " <<VectorToString(xc_puncturePoints[xc_i])<< std::endl;
allens's avatar
allens committed
977
            }
978 979

          if (seed == NULL)
allens's avatar
allens committed
980
            {
981 982
              seed = poincare_ic;
              poincare_ic->src_seed_ic = poincare_ic;
allens's avatar
allens committed
983
            }
984

985 986 987 988 989 990 991
          std::vector<avtPoincareIC *> *children = poincare_ic->src_rational_ic->properties.children;

          if ( !BracketIsValid( poincare_ic ) )
            {         
            if (2 <= RATIONAL_DEBUG)
              {
                std::cerr << "LINE " << __LINE__ << "  " << "Deleting ID: "<<_a->id << std::endl;
992
                std::cerr << "LINE " << __LINE__ << "  " << "Deleting ID: "<<_b->id << std::endl;
993 994 995 996 997 998 999 1000 1001
                std::cerr << "LINE " << __LINE__ << "  " << "Deleting ID: "<<_c->id << std::endl;
              }    

            if (_a->id != seed->id)
              ids_to_delete.push_back(_a->id);
            else
              _a->properties.searchState = FieldlineProperties::WAITING_SEED;
            if (_b->id != seed->id)
              ids_to_delete.push_back(_b->id);
allens's avatar
allens committed
1002 1003
            else
              _b->properties.searchState = FieldlineProperties::WAITING_SEED;