avtPoincareIC.h 13.2 KB
Newer Older
1 2
/*****************************************************************************
*
brugger's avatar
 
brugger committed
3
* Copyright (c) 2000 - 2013, Lawrence Livermore National Security, LLC
4 5 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 40 41 42 43 44 45 46
* Produced at the Lawrence Livermore National Laboratory
* LLNL-CODE-442911
* 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.
*
*****************************************************************************/

// ************************************************************************* //
//                            avtPoincareIC.h                              //
// ************************************************************************* //

#ifndef AVT_POINCARE_IC_H
#define AVT_POINCARE_IC_H

#include <avtStateRecorderIntegralCurve.h>
47 48
#include <map>

49 50
class vtkObject;

51 52 53 54 55 56 57 58 59 60
// ****************************************************************************
//  Class: avtPoincareIC
//
//  Purpose:
//      A derived type of avtStateRecorderIntegralCurve.  This class 
//      decides how to terminate a streamline.
//
//  Programmer: Hank Childs
//  Creation:   October 4, 2010
//
hrchilds's avatar
 
hrchilds committed
61 62 63 64 65
//  Modifications:
//
//    Hank Childs, Sun Dec  5 10:18:13 PST 2010
//    Add avtIVPField to CheckForTermination.
//
hrchilds's avatar
hrchilds committed
66 67 68 69
//   David Camp, Wed Mar  7 10:43:07 PST 2012
//   Added a Serialize flag to the arguments. This is to support the restore
//   ICs code.
//
70 71
// ****************************************************************************

72 73 74
#ifndef POINCARE_FIELDLINE_PROPERTIES_H
#define POINCARE_FIELDLINE_PROPERTIES_H

allens's avatar
allens committed
75 76
class avtPoincareIC;

77 78 79 80 81 82 83 84 85 86 87 88 89 90
class WindingPair {

public:
  WindingPair( unsigned int t, unsigned int p )
  {
    toroidal = t;
    poloidal = p;
  };

  unsigned int toroidal;
  unsigned int poloidal;
};


91 92 93 94 95 96 97 98
class FieldlineProperties {

public:

  FieldlineProperties()
  {
    type = FieldlineProperties::UNKNOWN_TYPE;

99 100 101
    analysisMethod = FieldlineProperties::DEFAULT_METHOD;
    searchState    = FieldlineProperties::NO_SEARCH;
    analysisState  = FieldlineProperties::UNKNOWN_ANALYSIS;
allens's avatar
allens committed
102

103 104 105 106
    source = FieldlineProperties::UNKNOWN_TYPE;
    
    iteration = 0;

allens's avatar
allens committed
107 108
    safetyFactor = 0;

109 110 111
    toroidalWinding = 0;
    poloidalWinding = 0;

allens's avatar
allens committed
112
    toroidalWindingP = 0;
allens's avatar
allens committed
113 114 115 116
    poloidalWindingP = 0;

    toroidalResonance = 0;
    poloidalResonance = 0;
117 118 119

    windingGroupOffset = 0;
    islands = 0;
allens's avatar
allens committed
120
    islandGroups = 0;
121 122 123

    nnodes  = 0;
    
allens's avatar
allens committed
124 125
    rationalSurfaceTolerance = 0;

126
    maxPunctures      = 0;
127
    numPunctures      = 0;
128
    nPuncturesNeeded  = 0;
allens's avatar
allens committed
129 130 131

    parentOPointIC = 0;
    childOPointIC = 0;
132 133 134 135
  };

enum FieldlineType { UNKNOWN_TYPE  = 0,

allens's avatar
allens committed
136 137 138 139 140
                     PERIODIC = 0x0010,
                     RATIONAL = 0x0010,

                     O_POINT  = 0x0011,
                     X_POINT  = 0x0012,
141 142

                     
allens's avatar
allens committed
143 144 145
                     QUASI_PERIODIC = 0x0060,
                     IRRATIONAL     = 0x0060,

146
                     FLUX_SURFACE   = 0x0020,
allens's avatar
allens committed
147 148

                     ISLAND_CHAIN                    = 0x0040,
allens's avatar
allens committed
149

allens's avatar
allens committed
150 151
                     ISLAND_PRIMARY_CHAIN            = 0x0040,
                     ISLAND_SECONDARY_CHAIN          = 0x0042,
allens's avatar
allens committed
152

allens's avatar
allens committed
153 154
                     ISLAND_PRIMARY_SECONDARY_AXIS   = 0x0041,
                     ISLAND_SECONDARY_SECONDARY_AXIS = 0x0043,
155 156 157
                     
                     CHAOTIC = 30 };
  
158 159
enum AnalysisMethod { UNKNOWN_METHOD = 0,
                      DEFAULT_METHOD = 1,
allens's avatar
allens committed
160

161
                      RATIONAL_SEARCH   = 10,
162
                      RATIONAL_MINIMIZE };
163 164

enum AnalysisState { UNKNOWN_ANALYSIS = 0,
165

allens's avatar
allens committed
166 167 168
                     COMPLETED  = 2,
                     TERMINATED = 3,
                     
allens's avatar
allens committed
169 170
                     OVERRIDE = 5,

allens's avatar
allens committed
171 172
                     DELETE   = 7,

173 174
                     ADDING_POINTS = 10,

allens's avatar
allens committed
175 176
                     ADD_O_POINTS = 15,
                     ADD_X_POINTS = 16,
177

178 179 180
                     ADD_BOUNDARY_POINT = 17 };

enum SearchState { UNKNOWN_SEARCH = 0,
181

182
                   NO_SEARCH = 1,
allens's avatar
allens committed
183

184 185 186
                   ////// Code for island width search
                   ISLAND_O_POINT,
                   ISLAND_BOUNDARY_SEARCH,
allens's avatar
allens committed
187

188
                   ////// Code for island width search
allens's avatar
allens committed
189
                      
190 191 192 193 194
                   ////// Code for rational surface search
                   ORIGINAL_RATIONAL = 100,
                   SEARCHING_SEED,
                   WAITING_SEED,
                   FINISHED_SEED,
195
                   DEAD_SEED,
196 197 198 199 200 201
                   MINIMIZING_A      = 105,  // Used to bracket the minimum
                   MINIMIZING_B,
                   MINIMIZING_C,
                   MINIMIZING_X0     = 110, // Used for Golden search routine
                   MINIMIZING_X1,
                   MINIMIZING_X2,
202
                   MINIMIZING_X3
allens's avatar
allens committed
203 204
                      ////// Code for rational surface search
};
allens's avatar
allens committed
205 206


207 208 209 210 211 212
public:

  FieldlineType type;

  FieldlineType source;

allens's avatar
allens committed
213 214 215
  ////// Code for rational surface search
  AnalysisMethod analysisMethod;

216 217
  AnalysisState analysisState;

allens's avatar
allens committed
218
  ////// Code for rational surface search
219
  SearchState searchState;
allens's avatar
allens committed
220

221 222
  unsigned int iteration;

allens's avatar
allens committed
223 224
  double safetyFactor;

allens's avatar
allens committed
225
  // Base number of transits
226 227
  unsigned int toroidalWinding;
  unsigned int poloidalWinding;
allens's avatar
allens committed
228

allens's avatar
allens committed
229
  // Secondary axis number of transits
allens's avatar
allens committed
230
  unsigned int toroidalWindingP;
allens's avatar
allens committed
231
  unsigned int poloidalWindingP;
232

allens's avatar
allens committed
233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248
  // Number of transists when an island to get back to the initial puncture.
  unsigned int toroidalPeriod;
  unsigned int poloidalPeriod;

  // Resonance periods 

  // When a surface resonances equal 1

  // When a primary island the resonances equal the base number of
  // windings and the toroidal resonance and the base toroidal winding both
  // equal the number of islands.

  // When secondary islands the toroial resonance is the total number
  // of islands and the toroidal resonance divided by the base
  // toroidal winding equals the number of island on each group.

allens's avatar
allens committed
249 250
  unsigned int toroidalResonance;
  unsigned int poloidalResonance;
allens's avatar
allens committed
251

252 253 254
  std::vector< WindingPair > windingPairs;

  std::map< WindingPair, unsigned int > topWindingPairs;
255 256 257

  unsigned int windingGroupOffset;
  unsigned int islands;
allens's avatar
allens committed
258
  unsigned int islandGroups;
259

allens's avatar
allens committed
260
  // If a surface it's overlap is found geometrically
allens's avatar
allens committed
261
  // If an island (primary or secondary) toroidalPeriod / toroidalResonance
262 263 264
  float nnodes;

  unsigned int maxPunctures;
265
  unsigned int numPunctures;
266 267
  unsigned int nPuncturesNeeded;

allens's avatar
allens committed
268 269 270
  // Rational Surface periodicity measures
  double rationalSurfaceTolerance;

allens's avatar
allens committed
271 272 273 274 275 276 277
  // Seeds for islands
  avtVector lastSeedPoint;
  std::vector< avtVector > seedPoints;

  ////// Code for island width search
  bool pastFirstSearchFailure;
  double islandWidth;
allens's avatar
allens committed
278
  double searchBaseDelta;
allens's avatar
allens committed
279 280 281 282 283 284 285 286 287 288 289
  double searchDelta;
  double searchIncrement;
  double searchMagnitude;
  avtVector searchNormal;
  avtPoincareIC* parentOPointIC;
  avtPoincareIC* childOPointIC;

  unsigned int baseToroidalWinding;
  unsigned int basePoloidalWinding;

  ////// Code for island width search
290

allens's avatar
allens committed
291 292 293 294 295 296
  ////// Code for rational surface search
  // The rational points bounding the location of the minimization action
  avtVector rationalPt1;
  avtVector rationalPt2;

  std::vector< avtPoincareIC *> *children;
297 298

  avtVector srcPt;
allens's avatar
allens committed
299
  ////// Code for rational surface search
300 301 302
};
#endif

303 304 305 306 307
class IVP_API avtPoincareIC : public avtStateRecorderIntegralCurve
{
public:
    avtPoincareIC(unsigned char mask, const avtIVPSolver* model, 
                  Direction dir, const double& t_start, 
308 309 310
                  const avtVector &p_start,
                  const avtVector &v_start,
                  int ID);
311 312 313 314 315 316 317

    void          SetIntersectionCriteria(vtkObject *obj, int);

    avtPoincareIC();
    virtual ~avtPoincareIC();

    virtual void  Serialize(MemStream::Mode mode, MemStream &buff, 
hrchilds's avatar
hrchilds committed
318
                            avtIVPSolver *solver, SerializeFlags serializeFlags);
319 320 321 322 323

  protected:
    avtPoincareIC( const avtPoincareIC& );
    avtPoincareIC& operator=( const avtPoincareIC& );
    
allens's avatar
allens committed
324
    bool         IntersectPlane(const avtVector &p0, const avtVector &p1);
325 326

  public:
hrchilds's avatar
 
hrchilds committed
327
    virtual bool CheckForTermination(avtIVPStep& step, avtIVPField *);
328

allens's avatar
allens committed
329 330 331
    bool         TerminatedBecauseOfMaxIntersections(void) 
                            { return terminatedBecauseOfMaxIntersections; };

332 333 334 335
    // Intersection points.
    bool   intersectionsSet;
    int    maxIntersections;
    int    numIntersections;
allens's avatar
allens committed
336
    double intersectPlaneEq[4]; // Typically the Y=0 plane i.e. 0, 1, 0
allens's avatar
allens committed
337
    bool   terminatedBecauseOfMaxIntersections;
338

allens's avatar
allens committed
339 340 341
    // These are the fieldline points as stripped out of the IC
    // proper.  They are stored here for convience so the analysis can
    // be done without schleping the whole integral curve around.
342 343
    std::vector<avtVector> points;

allens's avatar
allens committed
344
    // The fieldline properties as returned from the analysis library.
345
    FieldlineProperties properties;
346

allens's avatar
allens committed
347 348

    ////// Code for rational surface search
349 350
    avtPoincareIC *src_seed_ic;
    avtPoincareIC *src_rational_ic;
allens's avatar
allens committed
351 352 353 354 355 356 357 358 359 360 361

    // If this curve is minimizing, keep track of 'a' and 'c' (this is 'b')
    avtPoincareIC *a_IC;
    avtPoincareIC *b_IC;
    avtPoincareIC *c_IC;
    // Golden Search catches X0. X1, X2 and X3 must all have had integration done
    avtPoincareIC *GS_x1;
    avtPoincareIC *GS_x2;
    avtPoincareIC *GS_x3;
    ////// Code for rational surface search
};
allens's avatar
allens committed
362 363 364

// ostream operators for FieldlineProperties::AnalysisState's enum types
inline std::ostream& operator<<( std::ostream& out, 
365
                                 FieldlineProperties::AnalysisState state )
allens's avatar
allens committed
366
{
367
    switch( state )
allens's avatar
allens committed
368
    {
369 370
    case FieldlineProperties::UNKNOWN_ANALYSIS:
        return out << "UNKNOWN_ANALYSIS";
allens's avatar
allens committed
371 372 373 374
    case FieldlineProperties::OVERRIDE:
        return out << "OVERRIDE";
    case FieldlineProperties::ADDING_POINTS:
        return out << "ADDING_POINTS";
allens's avatar
allens committed
375 376
    case FieldlineProperties::ADD_O_POINTS:
        return out << "ADD_O_POINTS";
377 378
    case FieldlineProperties::ADD_BOUNDARY_POINT:
        return out << "ADD_BOUNDARY_POINT";
allens's avatar
allens committed
379 380 381 382 383
    case FieldlineProperties::COMPLETED:
        return out << "COMPLETED";
    case FieldlineProperties::TERMINATED:
        return out << "TERMINATED";
    default:
384
        return out << "???????";
allens's avatar
allens committed
385 386
    }
}
allens's avatar
allens committed
387

388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405
// ostream operators for FieldlineProperties::AnalysisState's enum types
inline std::ostream& operator<<( std::ostream& out, 
                                 FieldlineProperties::SearchState state )
{
    switch( state )
    {
    case FieldlineProperties::UNKNOWN_SEARCH:
        return out << "UNKNOWN_SEARCH";
    case FieldlineProperties::NO_SEARCH:
        return out << "NO_SEARCH";
    case FieldlineProperties::ISLAND_O_POINT:
        return out << "ISLAND_O_POINT";
    case FieldlineProperties::ISLAND_BOUNDARY_SEARCH:
        return out << "ISLAND_BOUNDARY_SEARCH";
    default:
        return out << "???????";
    }
}
allens's avatar
allens committed
406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432

inline std::ostream& operator<<( std::ostream& out, 
                                  FieldlineProperties::FieldlineType type )
{
    switch( type )
    {
    case FieldlineProperties::UNKNOWN_TYPE:
        return out << "UNKNOWN_TYPE";
    case FieldlineProperties::RATIONAL:
        return out << "RATIONAL";
    case FieldlineProperties::O_POINT:
        return out << "O_POINT";
    case FieldlineProperties::X_POINT:
        return out << "X_POINT";
    case FieldlineProperties::FLUX_SURFACE:
        return out << "FLUX_SURFACE";
    case FieldlineProperties::ISLAND_PRIMARY_CHAIN:
        return out << "ISLAND_PRIMARY_CHAIN";
    case FieldlineProperties::ISLAND_SECONDARY_CHAIN:
        return out << "ISLAND_SECONDARY_CHAIN";
    case FieldlineProperties::ISLAND_PRIMARY_SECONDARY_AXIS:
        return out << "ISLAND_PRIMARY_SECONDARY_AXIS";
    case FieldlineProperties::ISLAND_SECONDARY_SECONDARY_AXIS:
        return out << "ISLAND_SECONDARY_SECONDARY_AXIS";
    case FieldlineProperties::CHAOTIC:
        return out << "CHAOTIC";
    default:
433
        return out << "???????";
allens's avatar
allens committed
434 435
    }
}
436
#endif