avtM3DC1FileFormat.C 79.1 KB
Newer Older
1
2
3
// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.
4
5
6
7
8
9

// ************************************************************************* //
//                            avtM3DC1FileFormat.C                           //
// ************************************************************************* //

#include <avtM3DC1FileFormat.h>
Allen Sanderson's avatar
Allen Sanderson committed
10
#include <avtM3DC1Field.h>
11
12
13
14
15
16
17
18

#include <string>

#include <vtkIntArray.h>
#include <vtkFloatArray.h>
#include <vtkDoubleArray.h>
#include <vtkUnstructuredGrid.h>
#include <vtkTriangle.h>
Allen Sanderson's avatar
Allen Sanderson committed
19
#include <vtkWedge.h>
20
21
22
23

#include <avtDatabaseMetaData.h>

#include <DBOptionsAttributes.h>
24
25
#include <avtLogicalSelection.h>
//#include <avtSpatialBoxSelection.h>
26
27
#include <Expression.h>

28
29
#include <avtCallback.h>
#include <NonCompliantException.h>
30
#include <InvalidFilesException.h>
Allen Sanderson's avatar
Allen Sanderson committed
31
#include <InvalidVariableException.h>
32
33
#include <DebugStream.h>

34
#include <vtk_hdf5.h>
35
36
#include <visit-hdf5.h>

37

Allen Sanderson's avatar
Allen Sanderson committed
38
39
40
41
42
#define ELEMENT_SIZE_2D 7
#define SCALAR_SIZE_2D 20

#define ELEMENT_SIZE_3D 9
#define SCALAR_SIZE_3D 80
43
44
45
46
47
48
49
50
51
52
53
54
55


// ****************************************************************************
//  Method: avtM3DC1FileFormat constructor
//
//  Programmer: allen -- generated by xml2avt
//  Creation:   Fri Dec 4 15:04:15 PST 2009
//
// ****************************************************************************

avtM3DC1FileFormat::avtM3DC1FileFormat(const char *filename,
                                       DBOptionsAttributes *readOpts)
  : avtMTSDFileFormat(&filename, 1),
56
    processDataSelections(false), haveReadWholeData(true),
57
    m_filename(filename),
58
    m_refinement(2), m_dataLocation(AVT_NODECENT)
59
60
61
{
    if (readOpts != NULL) {
      for (int i=0; i<readOpts->GetNumberOfOptions(); ++i) {
Allen Sanderson's avatar
Allen Sanderson committed
62
        if (readOpts->GetName(i) == "Mesh refinement")
Allen Sanderson's avatar
Allen Sanderson committed
63
64
          m_refinement = readOpts->GetEnum("Mesh refinement");
        else if (readOpts->GetName(i) == "Linear mesh data location") 
65
        {
Allen Sanderson's avatar
Allen Sanderson committed
66
          int dataLocation = readOpts->GetEnum("Linear mesh data location");
67
68
69
70
71
72
          
          if( dataLocation == 0 )
            m_dataLocation = AVT_NODECENT;
          else if( dataLocation == 1 )
            m_dataLocation = AVT_ZONECENT;
        }
73
74
75
        else if (readOpts->GetName(i) == "Process Data Selections in the Reader")
          processDataSelections =
            readOpts->GetBool("Process Data Selections in the Reader");
76
77
78
      }
    }

Allen Sanderson's avatar
Allen Sanderson committed
79
    if( m_refinement < 0 )      m_refinement = 0;
Allen Sanderson's avatar
Allen Sanderson committed
80
    else if( m_refinement > 4 ) m_refinement = 4;
81

Allen Sanderson's avatar
Allen Sanderson committed
82
83
84
85
86
87
88
89
    if( processDataSelections && m_refinement )
    {
        EXCEPTION2( NonCompliantException, "M3DC1 Mesh",
                    "Can not subselect a refined mesh. " +
                    std::string(" Set the mesh refinement to 1 or disable") +
                    std::string(" processing data selections in the reader.") );
    }

90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
    LoadFile();
}


// ****************************************************************************
//  Method: avtEMSTDFileFormat::GetNTimesteps
//
//  Purpose:
//      Tells the rest of the code how many timesteps there are in this file.
//
//  Programmer: allen -- generated by xml2avt
//  Creation:   Fri Dec 4 15:04:15 PST 2009
//
// ****************************************************************************

int
avtM3DC1FileFormat::GetNTimesteps(void)
{
Kathleen Bonnell's avatar
Kathleen Bonnell committed
108
    return (int)m_times.size();
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
}


// ****************************************************************************
//  Method: avtM3DC1FileFormat::FreeUpResources
//
//  Purpose:
//      When VisIt is done focusing on a particular timestep, it asks that
//      timestep to free up any resources (memory, file descriptors) that
//      it has associated with it.  This method is the mechanism for doing
//      that.
//
//  Programmer: allen -- generated by xml2avt
//  Creation:   Fri Dec 4 15:04:15 PST 2009
//
// ****************************************************************************

void
avtM3DC1FileFormat::FreeUpResources(void)
{
    H5Fclose( m_fileID );
}


133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
// ****************************************************************************
//  Method: avtVsFileFormat::CanCacheVariable
//
//  Purpose:
//      To truly exercise the VS file format, we can't have VisIt caching
//      chunks of mesh and variables above the plugin.
//      
//  Programmer: Mark C. Miller 
//  Creation:   September 20, 2004 
//
// ****************************************************************************

bool
avtM3DC1FileFormat::CanCacheVariable(const char *var)
{
    // If processing the selections turn caching off.
    return !processDataSelections;

//    return haveReadWholeData;
}


// ****************************************************************************
//  Method: avtVsFileFormat::RegisterDataSelections
//
//  Purpose:
//      The Vs format can exploit some data selections so get them. 
//      
//  Programmer: Allen Sanderson
//  Creation:   March 4, 2011
//
// ****************************************************************************

void
avtM3DC1FileFormat::RegisterDataSelections(const std::vector<avtDataSelection_p> &sels,
                                        std::vector<bool> *selectionsApplied)
{
  selList     = sels;
  selsApplied = selectionsApplied;
}


// ****************************************************************************
//  Method: avtVsFileFormat::ProcessDataSelections
//
//  Purpose:
//      The Vs format can exploit some data selections so process them. 
//      
//  Programmer: Allen Sanderson
//  Creation:   March 4, 2011
//
// ****************************************************************************

bool
avtM3DC1FileFormat::ProcessDataSelections(int *mins, int *maxs, int *strides)
{
    bool retval = false;

    if( !processDataSelections )
      return retval;

    avtLogicalSelection composedSel;

196
    for (size_t i = 0; i < selList.size(); i++)
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
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
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
    {
        if (std::string(selList[i]->GetType()) == "Logical Data Selection")
        {
            avtLogicalSelection *sel = (avtLogicalSelection *) *(selList[i]);

            // overrwrite method-scope arrays with the new indexing
            composedSel.Compose(*sel);
            (*selsApplied)[i] = true;
            retval = true;
        }

        // Cannot handle avtSpatialBoxSelection without knowing the mesh.

//         else if (std::string(selList[i]->GetType()) == "Spatial Box Data Selection")
//         {
//             avtSpatialBoxSelection *sel =
//               (avtSpatialBoxSelection *) *(Sellist[i]);

//             double dmins[3], dmaxs[3];
//             sel->GetMins(dmins);
//             sel->GetMaxs(dmaxs);
//             avtSpatialBoxSelection::InclusionMode imode =
//                 sel->GetInclusionMode();

//             // we won't handle clipping of zones here
//             if ((imode != avtSpatialBoxSelection::Whole) &&
//                 (imode != avtSpatialBoxSelection::Partial))
//             {
//                 (*selsApplied)[i] = false;
//                 continue;
//             }

//             int imins[3], imaxs[3];
//             for (int j = 0; j < 3; j++)
//             {
//                 int imin = (int) dmins[j];
//                 if (((double) imin < dmins[j]) &&
//                     (imode == avtSpatialBoxSelection::Whole))
//                     imin++;
                
//                 int imax = (int) dmaxs[j];
//                 if (((double) imax < dmaxs[j]) &&
//                     (imode == avtSpatialBoxSelection::Partial))
//                     imax++;

//                 imins[j] = imin;
//                 imaxs[j] = imax;
//             }

//             avtLogicalSelection newSel;
//             newSel.SetStarts(imins);
//             newSel.SetStops(imaxs);

//             composedSel.Compose(newSel);
//             (*selsApplied)[i] = true;
//             retval = true;
//         }
        else
        {
            // indicate we won't handle this selection
            (*selsApplied)[i] = false;
        }
    }

    composedSel.GetStarts(mins);
    composedSel.GetStops(maxs);
    composedSel.GetStrides(strides);

265
    // If the user is a dumb ass and selects a dimension lower than
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
    // the actual dimension the min, max, and stride will be zero. So
    // fix it to be the full bounds and a stride of 1.
    for (int i = 0; i < 3; i++)
    {
      if( strides[i] == 0 )
      {
        mins[i] = 0;
        maxs[i] = -1;
        strides[i] = 1;
      }
    }

    return retval;
}


282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
// ****************************************************************************
//  Method: avtM3DC1FileFormat::PopulateDatabaseMetaData
//
//  Purpose:
//      This database meta-data object is like a table of contents for the
//      file.  By populating it, you are telling the rest of VisIt what
//      information it can request from you.
//
//  Programmer: allen -- generated by xml2avt
//  Creation:   Fri Dec 4 15:04:15 PST 2009
//
// ****************************************************************************

void
avtM3DC1FileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *md,
                                             int timeState)
{
    avtMeshMetaData *mmd;
300
    avtMeshType meshType = AVT_UNSTRUCTURED_MESH;
301
    int nblocks = 1;
302
303
304
305
306
    int block_origin = 0;
    int cell_origin = 0;
    int group_origin = 0;
    int spatial_dimension = 3;
    int topological_dimension = 3;
307
    //double *extents = NULL;
308
    int bounds[3] = {nelms, 0, 0};
309
310
311

    char level[4];

Allen Sanderson's avatar
Allen Sanderson committed
312
    if( m_refinement != 0 )
313
    {
Allen Sanderson's avatar
Allen Sanderson committed
314
315
      sprintf( level, "_%d", m_refinement );
    }
316
    else
317
      level[0] = '\0';
318
319
320


    // Original meshes for the user to see.
321
322
323
324
325
    mmd = new avtMeshMetaData("equilibrium/mesh",
                              nblocks, block_origin,
                              cell_origin, group_origin,
                              spatial_dimension, topological_dimension,
                              meshType);
Allen Sanderson's avatar
Allen Sanderson committed
326
    
327
328
329
330
331
332
333
334
335
336
337
338
339
340
    mmd->SetBounds( bounds );
    mmd->SetNumberCells( nelms );
    md->Add(mmd);

    mmd = new avtMeshMetaData("mesh",
                              nblocks, block_origin,
                              cell_origin, group_origin,
                              spatial_dimension, topological_dimension,
                              meshType);

    mmd->SetBounds( bounds );
    mmd->SetNumberCells( nelms );
    md->Add(mmd);

341
342
    // Populate the scalar field vars that will be interpolate onto a
    // refined mesh.
343
    for ( size_t i = 0; i < m_fieldVarNames.size(); ++i )
344
    {
345
346
      std::string varname = "equilibrium/" + m_fieldVarNames[i];
      std::string meshname = std::string("equilibrium/mesh") + std::string(level);
Allen Sanderson's avatar
Allen Sanderson committed
347
      
348
      AddScalarVarToMetaData( md, varname, meshname, m_dataLocation );
Allen Sanderson's avatar
Allen Sanderson committed
349
      
350
      meshname = std::string("mesh") + std::string(level);
351
352
353
      AddScalarVarToMetaData( md, m_fieldVarNames[i], meshname, m_dataLocation );
    }

354
355
356
    // For now the mesh is the same mesh as the original mesh because
    // of needing it for the integration.
    AddVectorVarToMetaData( md, "B_C1_Elements",
357
                            std::string("mesh"),
358
359
360
361
                            AVT_ZONECENT, 3);
    
    // Interpolated on to a mesh for visualization only.
    AddVectorVarToMetaData( md, "B_Interpolated",
362
                            std::string("mesh") + std::string(level),
363
364
                            m_dataLocation, 3 );
    
365
    // Hidden refined meshes for working with the interpolated data
Allen Sanderson's avatar
Allen Sanderson committed
366
    if( m_refinement )
367
    {
Allen Sanderson's avatar
Allen Sanderson committed
368
369
      int nLevels = m_refinement + 1;

Allen Sanderson's avatar
Allen Sanderson committed
370
      nblocks = nelms * nLevels * nLevels;
371
372

      mmd =
373
        new avtMeshMetaData(std::string("equilibrium/mesh") + std::string(level),
374
375
                            nblocks, block_origin,
                            cell_origin, group_origin,
376
                            spatial_dimension, topological_dimension, meshType);
377
378
379
380
      mmd->hideFromGUI = true;
      md->Add(mmd);

      mmd =
381
        new avtMeshMetaData(std::string("mesh") + std::string(level),
382
383
                            nblocks, block_origin,
                            cell_origin, group_origin,
384
                            spatial_dimension, topological_dimension, meshType);
385
386
387
388
389
390
391
392
393
394
      mmd->hideFromGUI = true;
      md->Add(mmd);
    }

    nblocks = nelms;    

    // Hidden meshes for working with the elements directly
    mmd =
      new avtMeshMetaData("hidden/equilibrium/mesh",
                          nblocks, block_origin, cell_origin, group_origin,
395
                          spatial_dimension, topological_dimension, meshType);
396
397
398
399
400
401
    mmd->hideFromGUI = true;
    md->Add(mmd);

    mmd =
      new avtMeshMetaData("hidden/mesh",
                          nblocks, block_origin, cell_origin, group_origin,
402
                          spatial_dimension, topological_dimension, meshType);
403
404
405
406
407
    mmd->hideFromGUI = true;
    md->Add(mmd);


    // Hidden scalar header vars.
408
    for ( size_t i = 0; i < m_scalarVarNames.size(); ++i )
409
410
411
412
    {
      avtScalarMetaData *smd =
        new avtScalarMetaData("hidden/" + m_scalarVarNames[i],
                              "hidden/equilibrium/mesh", AVT_ZONECENT);
Allen Sanderson's avatar
Allen Sanderson committed
413
      
414
415
416
      smd->hideFromGUI = true;
      md->Add(smd);
    }
Allen Sanderson's avatar
Allen Sanderson committed
417
    
418
419
420
    // Add the elements so we have access to them for the interpolation
    avtVectorMetaData *amd =
      new avtVectorMetaData("hidden/equilibrium/elements",
Allen Sanderson's avatar
Allen Sanderson committed
421
422
423
                            "hidden/equilibrium/mesh",
                            AVT_ZONECENT, element_size);
    
424
425
    amd->hideFromGUI = true;
    md->Add(amd);
Allen Sanderson's avatar
Allen Sanderson committed
426
    
Allen Sanderson's avatar
bug fix    
Allen Sanderson committed
427
    amd =
Allen Sanderson's avatar
Allen Sanderson committed
428
429
      new avtVectorMetaData("hidden/elements", "hidden/mesh",
                            AVT_ZONECENT, element_size);
430
431
432
433
434

    amd->hideFromGUI = true;
    md->Add(amd);

    // Hidden array field vars so we have access to them for the interpolation
435
    std::string varname;
Allen Sanderson's avatar
Allen Sanderson committed
436

437
    for ( size_t i = 0; i < m_fieldVarNames.size(); ++i )
438
    {
Allen Sanderson's avatar
Allen Sanderson committed
439
      varname = "hidden/equilibrium/" + m_fieldVarNames[i];
440
      amd = new avtVectorMetaData(varname, "hidden/equilibrium/mesh",
Allen Sanderson's avatar
Allen Sanderson committed
441
442
                                  AVT_ZONECENT, scalar_size);
      
443
444
      amd->hideFromGUI = true;
      md->Add(amd);
Allen Sanderson's avatar
Allen Sanderson committed
445
    
446
447
      varname = "hidden/" + m_fieldVarNames[i];
      amd = new avtVectorMetaData(varname, "hidden/mesh",
Allen Sanderson's avatar
Allen Sanderson committed
448
                                 AVT_ZONECENT, scalar_size);
449
450
451
452

      amd->hideFromGUI = true;
      md->Add(amd);
    }
Allen Sanderson's avatar
Allen Sanderson committed
453
454
455
456
457
458

    md->SetCyclesAreAccurate(true);
    md->SetCycles( m_cycles );

    md->SetTimesAreAccurate(true);
    md->SetTimes( m_times );
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
}


// ****************************************************************************
//  Method: avtM3DC1FileFormat::GetElements
//
//  Purpose:
//      Gets the elements assoicated for the C1 mesh.
//
//  Arguments:
//      timestate   The index of the timestate.  If GetNTimesteps returned
//                  'N' time steps, this is guaranteed to be between 0 and N-1.
//      meshname    The name of the mesh of interest.  This can be ignored if
//                  there is only one mesh.
//
//  Programmer: allen -- generated by xml2avt
//  Creation:   Fri Dec 4 15:04:15 PST 2009
//
// ****************************************************************************

float *
avtM3DC1FileFormat::GetElements(int timestate, const char *meshname)
{
  char meshStr[64];

  // Parse the mesh naem into the hdf5 group name.
  if( strncmp(meshname, "equilibrium/mesh", 16 ) == 0 )
  {
    sprintf( meshStr, "/equilibrium/mesh" );
  } else if( strncmp(meshname, "mesh", 4 ) == 0 ) {
    sprintf( meshStr, "/time_%03d/mesh", timestate );
  } else
491
    EXCEPTION2( NonCompliantException, "M3DC1 Element Name Lookup",
492
                "Element '" + std::string(meshStr) + "' was not found." );
493
494
495
496

  // Open the group.
  hid_t meshId = H5Gopen( m_fileID, meshStr, H5P_DEFAULT);
  if ( meshId < 0 )
497
    EXCEPTION2( NonCompliantException, "M3DC1 Group Open",
498
                "Group '" + std::string(meshStr) + "' was not found." );
499
500
501
502
  
  // Read in the mesh information.
  int nElements;
  if ( ! ReadAttribute( meshId, "nelms", &nElements ) )
503
504
    EXCEPTION2( NonCompliantException, "M3DC1 Attribute Reader",
                "Attribute 'nelms' was not found or was the wrong type." );
505
506
  
  if( nElements != nelms )
507
508
    EXCEPTION2( NonCompliantException, "M3DC1 Element Check",
                "Time step 'nelms' does not match equilibrium 'nelms'" );
509
510
511
512
513
 
  // Open the dataset and space info for the elements.  
  hid_t datasetId = H5Dopen(meshId, "elements", H5P_DEFAULT);
  hid_t spaceId = H5Dget_space(datasetId);
  size_t rank = H5Sget_simple_extent_ndims(spaceId);
514
515
  std::vector<hsize_t> sdim(rank);
  H5Sget_simple_extent_dims(spaceId, &sdim[0], NULL);
516
517

  // Sanity check.  
518
  if( rank != 2 || (size_t)sdim[0] != nelms || sdim[1] != element_size )
519
  {
520
521
      EXCEPTION2( NonCompliantException, "M3DC1 Element Check",
                  "The number of elements or the element size does not match" );          
522
523
524
525
526
  }

  // Memory for the elements.  
  float *elements = new float[sdim[0]*sdim[1]];
  if( elements == 0 )
527
528
    EXCEPTION2( NonCompliantException, "M3DC1 Memory Allocation",
                "CAN NOT ALLOCATE MEMORY" );
529
530
531
532
533

  // Read in the elements.
  H5Dread( datasetId,
           H5T_NATIVE_FLOAT, H5S_ALL, spaceId, H5P_DEFAULT, elements );

534
  H5Sclose(spaceId);
535
536
537
538
539
540
541
542
543
544
545
546
547
548
  H5Dclose(datasetId);
  H5Gclose( meshId );

  return elements;
}

// ****************************************************************************
//  Method: avtM3DC1FileFormat::GetMeshPoints
//
//  Purpose: Gets the mesh points associated with this file.  The mesh
//      is returned as vtkPoints.
//
//  Arguments:
//      elements        The original C1 elements 
Allen Sanderson's avatar
Allen Sanderson committed
549
//      refinement      The amount of mesh refinement.
550
551
552
553
554
555
556
//
//  Programmer: allen -- generated by xml2avt
//  Creation:   Fri Dec 4 15:04:15 PST 2009
//
// ****************************************************************************

vtkPoints *
Allen Sanderson's avatar
Allen Sanderson committed
557
avtM3DC1FileFormat::GetMeshPoints( float *elements, int refinement )
558
{
Allen Sanderson's avatar
Allen Sanderson committed
559
560
561
  // Inital number of triangles before refinement.
  int npts = nelms * nvertices;
  int ncoords = nvertices * 3;
562
563
564

  // VTK structure for holding the mesh points. 
  vtkPoints *vtkPts = vtkPoints::New();
Allen Sanderson's avatar
Allen Sanderson committed
565
  vtkPts->SetNumberOfPoints( npts );
566
567
  float *pts = (float *) vtkPts->GetVoidPointer(0);

Allen Sanderson's avatar
Allen Sanderson committed
568
569
570
571
572
573
574
  float phi = 0;

  // Pointer for fast indexing through the elements.
  float *element = elements;

//   for( int i=0; i<nelms/2; ++i )
//     element += element_size;
575

Allen Sanderson's avatar
Allen Sanderson committed
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
  for( int i=0; i<nelms; ++i )
  {
    // The element values from Jardin, JCP 200:133 (2004)
    float a     = element[0];
    float b     = element[1];
    float c     = element[2];
    float theta = element[3];
    float x     = element[4];
    float z     = element[5];

    // The seventh value (call it "e") specifies which edges of the
    // element lie on the domain boundary.  The first edge is the
    // edge between the first and second nodes, the second edge
    // joins nodes 2 and 3, and the third edge joins nodes 3 and 1.
    // If edge n lies on the domain boundary, than the nth bit of e
    // is set.  For example, e=0 means no edges are on the domain
    // boundary, e=5 means the first and third edges are on the
    // domain boundary.
594
    //unsigned int e = element[6];
Allen Sanderson's avatar
Allen Sanderson committed
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610

    if( element_dimension == 3 )
      phi = element[8];

    // The three points in the phi plane that define the triangle mesh element.
    pts[0] = x;
    pts[1] = phi;
    pts[2] = z;

    pts[3] = x + (a+b)*cos(theta);
    pts[4] = phi;
    pts[5] = z + (a+b)*sin(theta);
    
    pts[6] = x + b*cos(theta) - c*sin(theta);
    pts[7] = phi;
    pts[8] = z + b*sin(theta) + c*cos(theta);
611

Allen Sanderson's avatar
Allen Sanderson committed
612
613
614
    // The three points in the phi+d plane that define the second part of
    // the wedge element.
    if( element_dimension == 3 )
615
    {
Allen Sanderson's avatar
Allen Sanderson committed
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
      float d = element[7];

      pts[ 9] = x;
      pts[10] = phi + d;
      pts[11] = z;

      pts[12] = x + (a+b)*cos(theta);
      pts[13] = phi + d;
      pts[14] = z + (a+b)*sin(theta);
    
      pts[15] = x + b*cos(theta) - c*sin(theta);
      pts[16] = phi + d;
      pts[17] = z + b*sin(theta) + c*cos(theta);
    }

    // Increment for fast indexing.
    element += element_size;
    pts += ncoords;
  } 
Allen Sanderson's avatar
Allen Sanderson committed
635
636
637
638

  // Now do any needed refinement. Currently the flag is 0-4. The user
  // sees 1-5. With 1 being no refinement aka linear.
  if( refinement )
639
  {
Allen Sanderson's avatar
Allen Sanderson committed
640
641
    int nLevels = refinement + 1;

642
643
644
645
646
    // Pointer to the current mesh.
    float *pts = (float *) vtkPts->GetVoidPointer(0);

    // VTK structure for holding the refined mesh points. 
    vtkPoints *rf_vtkPts = vtkPoints::New();
Allen Sanderson's avatar
Allen Sanderson committed
647
648
649
650
651

    rf_vtkPts->SetNumberOfPoints( npts *
                                  (int) pow( (float) nLevels,
                                             (float) element_dimension) );

652
653
    float *rf_pts = (float *) rf_vtkPts->GetVoidPointer(0);

Allen Sanderson's avatar
Allen Sanderson committed
654
    if( element_dimension == 2 )
655
    {
Allen Sanderson's avatar
Allen Sanderson committed
656
      for( int i=0; i<npts; i+=nvertices )
657
      {
Allen Sanderson's avatar
Allen Sanderson committed
658
659
660
661
        // Refine the triangle based on barrycentric coordinates.
        float *A = &(pts[0]);
        float *B = &(pts[3]);
        float *C = &(pts[6]);
Allen Sanderson's avatar
Allen Sanderson committed
662

Allen Sanderson's avatar
Allen Sanderson committed
663
        for( int j=0; j<nLevels; ++j )
Allen Sanderson's avatar
Allen Sanderson committed
664
        {
Allen Sanderson's avatar
Allen Sanderson committed
665
666
          float a  = (float)  j      / (float) nLevels;
          float a1 = (float) (j + 1) / (float) nLevels;
Allen Sanderson's avatar
Allen Sanderson committed
667

Allen Sanderson's avatar
Allen Sanderson committed
668
669
670
671
          for( int k=0; k<nLevels-j; ++k )
          {
            float b   = (float) (nLevels-j-k)   / (float) nLevels;
            float b_1 = (float) (nLevels-j-k-1) / (float) nLevels;
Allen Sanderson's avatar
Allen Sanderson committed
672

Allen Sanderson's avatar
Allen Sanderson committed
673
674
            float c  = (float)  k    / (float) nLevels;
            float c1 = (float) (k+1) / (float) nLevels;
Allen Sanderson's avatar
Allen Sanderson committed
675

Allen Sanderson's avatar
Allen Sanderson committed
676
677
678
            rf_pts[0] = a * A[0] + b * B[0] + c * C[0];
            rf_pts[1] = a * A[1] + b * B[1] + c * C[1];
            rf_pts[2] = a * A[2] + b * B[2] + c * C[2];
Allen Sanderson's avatar
Allen Sanderson committed
679

Allen Sanderson's avatar
Allen Sanderson committed
680
681
682
            rf_pts[3] = a * A[0] + b_1 * B[0] + c1 * C[0];
            rf_pts[4] = a * A[1] + b_1 * B[1] + c1 * C[1];
            rf_pts[5] = a * A[2] + b_1 * B[2] + c1 * C[2];
Allen Sanderson's avatar
Allen Sanderson committed
683

Allen Sanderson's avatar
Allen Sanderson committed
684
685
686
            rf_pts[6] = a1 * A[0] + b_1 * B[0] + c * C[0];
            rf_pts[7] = a1 * A[1] + b_1 * B[1] + c * C[1];
            rf_pts[8] = a1 * A[2] + b_1 * B[2] + c * C[2];
Allen Sanderson's avatar
Allen Sanderson committed
687

Allen Sanderson's avatar
Allen Sanderson committed
688
            rf_pts += ncoords;
Allen Sanderson's avatar
Allen Sanderson committed
689

Allen Sanderson's avatar
Allen Sanderson committed
690
691
692
693
694
695
696
            if( k )
            {
              float c_1 = (float) (k-1) / (float) nLevels;

              rf_pts[0] = a * A[0] + b * B[0] + c * C[0];
              rf_pts[1] = a * A[1] + b * B[1] + c * C[1];
              rf_pts[2] = a * A[2] + b * B[2] + c * C[2];
Allen Sanderson's avatar
Allen Sanderson committed
697
            
Allen Sanderson's avatar
Allen Sanderson committed
698
699
700
              rf_pts[3] = a1 * A[0] + b_1 * B[0] + c * C[0];
              rf_pts[4] = a1 * A[1] + b_1 * B[1] + c * C[1];
              rf_pts[5] = a1 * A[2] + b_1 * B[2] + c * C[2];
Allen Sanderson's avatar
Allen Sanderson committed
701
            
Allen Sanderson's avatar
Allen Sanderson committed
702
703
704
              rf_pts[6] = a1 * A[0] + b * B[0] + c_1 * C[0];
              rf_pts[7] = a1 * A[1] + b * B[1] + c_1 * C[1];
              rf_pts[8] = a1 * A[2] + b * B[2] + c_1 * C[2];
Allen Sanderson's avatar
Allen Sanderson committed
705
            
Allen Sanderson's avatar
Allen Sanderson committed
706
707
              rf_pts += ncoords;
            }
Allen Sanderson's avatar
Allen Sanderson committed
708
709
710
          }
        }

Allen Sanderson's avatar
Allen Sanderson committed
711
        pts += ncoords;
712
      }
Allen Sanderson's avatar
Allen Sanderson committed
713
714
715
716
717
718
719
720
721
722
723
724
725
726
    }

    else //if( element_dimension == 3 )
    {
      for( int i=0; i<npts; i+=nvertices )
      {
        // Refine the triangle based on barrycentric coordinates.
        float *A = &(pts[0]);
        float *B = &(pts[3]);
        float *C = &(pts[6]);

        float *D = &(pts[9]);
        float *E = &(pts[12]);
        float *F = &(pts[15]);
727

Allen Sanderson's avatar
Allen Sanderson committed
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
        for( int j=0; j<nLevels; ++j )
        {
          float a  = (float)  j      / (float) nLevels;
          float a1 = (float) (j + 1) / (float) nLevels;
          
          for( int k=0; k<nLevels-j; ++k )
          {
            float b   = (float) (nLevels-j-k)   / (float) nLevels;
            float b_1 = (float) (nLevels-j-k-1) / (float) nLevels;
            
            float c  = (float)  k    / (float) nLevels;
            float c1 = (float) (k+1) / (float) nLevels;

            double left_pts[9], right_pts[9];

            left_pts[0] = a * A[0] + b * B[0] + c * C[0];
            left_pts[1] = a * A[1] + b * B[1] + c * C[1];
            left_pts[2] = a * A[2] + b * B[2] + c * C[2];

            left_pts[3] = a * A[0] + b_1 * B[0] + c1 * C[0];
            left_pts[4] = a * A[1] + b_1 * B[1] + c1 * C[1];
            left_pts[5] = a * A[2] + b_1 * B[2] + c1 * C[2];

            left_pts[6] = a1 * A[0] + b_1 * B[0] + c * C[0];
            left_pts[7] = a1 * A[1] + b_1 * B[1] + c * C[1];
            left_pts[8] = a1 * A[2] + b_1 * B[2] + c * C[2];


            right_pts[0] = a * D[0] + b * E[0] + c * F[0];
            right_pts[1] = a * D[1] + b * E[1] + c * F[1];
            right_pts[2] = a * D[2] + b * E[2] + c * F[2];

            right_pts[3] = a * D[0] + b_1 * E[0] + c1 * F[0];
            right_pts[4] = a * D[1] + b_1 * E[1] + c1 * F[1];
            right_pts[5] = a * D[2] + b_1 * E[2] + c1 * F[2];

            right_pts[6] = a1 * D[0] + b_1 * E[0] + c * F[0];
            right_pts[7] = a1 * D[1] + b_1 * E[1] + c * F[1];
            right_pts[8] = a1 * D[2] + b_1 * E[2] + c * F[2];

            for( int l=0; l<nLevels; ++l )
            {
              for( int ii=0; ii<9; ++ii )
                rf_pts[ii] = left_pts[ii] + (double) l *
                  (right_pts[ii] - left_pts[ii]) / (double) nLevels;

              rf_pts += ncoords/2;

              for( int ii=0; ii<9; ++ii )
                rf_pts[ii] = left_pts[ii] + (double) (l+1) *
                  (right_pts[ii] - left_pts[ii]) / (double) nLevels;

              rf_pts += ncoords/2;
            }

            if( k )
            {
              float c_1 = (float) (k-1) / (float) nLevels;

              left_pts[0] = a * A[0] + b * B[0] + c * C[0];
              left_pts[1] = a * A[1] + b * B[1] + c * C[1];
              left_pts[2] = a * A[2] + b * B[2] + c * C[2];
              
              left_pts[3] = a1 * A[0] + b_1 * B[0] + c * C[0];
              left_pts[4] = a1 * A[1] + b_1 * B[1] + c * C[1];
              left_pts[5] = a1 * A[2] + b_1 * B[2] + c * C[2];
              
              left_pts[6] = a1 * A[0] + b * B[0] + c_1 * C[0];
              left_pts[7] = a1 * A[1] + b * B[1] + c_1 * C[1];
              left_pts[8] = a1 * A[2] + b * B[2] + c_1 * C[2];
              
              right_pts[0] = a * D[0] + b * E[0] + c * F[0];
              right_pts[1] = a * D[1] + b * E[1] + c * F[1];
              right_pts[2] = a * D[2] + b * E[2] + c * F[2];
              
              right_pts[3] = a1 * D[0] + b_1 * E[0] + c * F[0];
              right_pts[4] = a1 * D[1] + b_1 * E[1] + c * F[1];
              right_pts[5] = a1 * D[2] + b_1 * E[2] + c * F[2];
              
              right_pts[6] = a1 * D[0] + b * E[0] + c_1 * F[0];
              right_pts[7] = a1 * D[1] + b * E[1] + c_1 * F[1];
              right_pts[8] = a1 * D[2] + b * E[2] + c_1 * F[2];

              for( int l=0; l<nLevels; ++l )
              {
                for( int ii=0; ii<9; ++ii )
                  rf_pts[ii] = left_pts[ii] + (double) l *
                    (right_pts[ii] - left_pts[ii]) / (double) nLevels;
                
                rf_pts += ncoords/2;
                
                for( int ii=0; ii<9; ++ii )
                  rf_pts[ii] = left_pts[ii] + (double) (l+1) *
                    (right_pts[ii] - left_pts[ii]) / (double) nLevels;
                
                rf_pts += ncoords/2;
              }
            }
          }
        }

        pts += ncoords;
      }
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
    }

    // Delete the old points.
    vtkPts->Delete();

    // Update the pointers to the new mesh
    vtkPts = rf_vtkPts;
  }

  return vtkPts;  
}



// ****************************************************************************
//  Method: avtM3DC1FileFormat::GetMesh
//
//  Purpose:
//      Gets the mesh associated with this file.  The mesh is returned as a
//      derived type of vtkDataSet (ie vtkRectilinearGrid, vtkStructuredGrid,
//      vtkUnstructuredGrid, etc).
//
//  Arguments:
//      timestate   The index of the timestate.  If GetNTimesteps returned
//                  'N' time steps, this is guaranteed to be between 0 and N-1.
//      meshname    The name of the mesh of interest.  This can be ignored if
//                  there is only one mesh.
//
//  Programmer: allen -- generated by xml2avt
//  Creation:   Fri Dec 4 15:04:15 PST 2009
//
// ****************************************************************************

vtkDataSet *
avtM3DC1FileFormat::GetMesh(int timestate, const char *meshname)
{
  const char *meshnamePtr = meshname;  
  char meshStr[64];

870
871
872
  bool haveDataSelections;
  int mins[3], maxs[3], strides[3];
  
Allen Sanderson's avatar
Allen Sanderson committed
873
874
875
876
877
878
879
880
881
882
883
884
885
  // Parse the mesh variable name to get the true mesh name and the
  // refiement level in the meshnamePtr.
  if( strncmp(meshnamePtr, "equilibrium/mesh", 16 ) == 0 )
  {
    meshnamePtr = &(meshnamePtr[16]);
    sprintf( meshStr, "equilibrium/mesh" );    
  }
  else if( strncmp(meshnamePtr, "mesh", 4 ) == 0 ) {
    meshnamePtr = &(meshnamePtr[4]);
    sprintf( meshStr, "mesh" );    
  }
  else
    EXCEPTION2( NonCompliantException, "M3DC1 Mesh Name",
886
                "Can not find '" + std::string(meshnamePtr) + "'" );
Allen Sanderson's avatar
Allen Sanderson committed
887

888
889
  bool refinedMesh = strlen(meshnamePtr);
  
Allen Sanderson's avatar
Allen Sanderson committed
890
891
892
  // Parse the mesh variable name to get the refinement level.
  int refinement;

893
  if( refinedMesh && atoi(&(meshnamePtr[1])) == m_refinement )
Allen Sanderson's avatar
Allen Sanderson committed
894
895
896
897
      refinement = m_refinement;
  else
    refinement = 0;

898
899
900
901
902
  // If hidden only one poloidal plane as the mesh will be used for
  // interpolation. There should also be no refinement.
  if( strncmp(meshname, "hidden/", 7 ) == 0 )
  {
    meshnamePtr = &(meshname[7]);
903
904

    mins[0] = 0;
905
    maxs[0] = nelms-1;
906
907
908
    strides[0] = 1;
      
    haveReadWholeData = true;
909
910
911
912
  }
  else
  {
    meshnamePtr = meshname;
913
914

    // Adjust for the data selections which are NODAL.
915
    if( (haveDataSelections = ProcessDataSelections(mins, maxs, strides) ))  /// TODO: check on assignment
916
    {
Allen Sanderson's avatar
Allen Sanderson committed
917
918
919
920
921
922
923
924
      if( refinement == 0 )
      {
        EXCEPTION2( NonCompliantException, "M3DC1 Mesh Name",
                    "Can not subselect a refined mesh '" +
                    std::string(meshnamePtr) +
                    std::string("'. Set the mesh refinement to 1.") );
      }

925
926
927
928
929
930
931
932
933
934
935
936
      debug5
        << "Have a logical cell selection for mesh  " << meshname << "  "
        << "(" << mins[0] << "," << maxs[0] << " stride " << strides[0] << ") "
        //      << "(" << mins[1] << "," << maxs[1] << " stride " << strides[1] << ") "
        //      << "(" << mins[2] << "," << maxs[2] << " stride " << strides[2] << ") "
        << std::endl;
      
      haveReadWholeData = false;
    }
    else
    {
      mins[0] = 0;
Allen Sanderson's avatar
Allen Sanderson committed
937
938

      if ( refinement == 0 )
Kathleen Bonnell's avatar
Kathleen Bonnell committed
939
        maxs[0] = nelms-1;
Allen Sanderson's avatar
Allen Sanderson committed
940
941
942
      else if( element_dimension == 2 )
        maxs[0] = nelms*pow((double) refinement+1.0, 2.0)-1;
      else if( element_dimension == 3 )
Allen Sanderson's avatar
Allen Sanderson committed
943
944
        maxs[0] = nelms*pow((double) refinement+1.0, 3.0)-1;

945
946
947
948
      strides[0] = 1;
      
      haveReadWholeData = true;
    }
949
950
951
952
953
  }
  
  // Get the C1 elements.
  float* elements = GetElements(timestate, meshStr);

954
  // Create a VTK grid for the mesh.
955
956
  vtkUnstructuredGrid *grid = vtkUnstructuredGrid::New();

Allen Sanderson's avatar
Allen Sanderson committed
957
  // Now get the points on the mesh based on the refinement level.
958
  vtkPoints *vtkPts = GetMeshPoints( elements, refinement );
959
960

  // Add the points to the VTK grid.
961
  //int npts = vtkPts->GetNumberOfPoints();
962
963
964
965
966
967
968
  grid->SetPoints( vtkPts );

  delete [] elements;

  // Add in the VTK cells. The connectivity is the same for all
  // triangles because each triangle is defined by three points that
  // are not unique. 
969

970
971
  grid->Allocate((maxs[0]-mins[0])/strides[0]+1);

Allen Sanderson's avatar
Allen Sanderson committed
972
  if( element_dimension == 2 )
973
  {
Allen Sanderson's avatar
Allen Sanderson committed
974
    vtkTriangle *tri = vtkTriangle::New();
975
976

    for( int i=mins[0]; i<=maxs[0]; i+=strides[0] )
Allen Sanderson's avatar
Allen Sanderson committed
977
    {
978
979
980
981
      unsigned int index = i * nvertices;
      tri->GetPointIds()->SetId( 0, index   );
      tri->GetPointIds()->SetId( 1, index+1 );
      tri->GetPointIds()->SetId( 2, index+2 );
Allen Sanderson's avatar
Allen Sanderson committed
982
983
984
      
      grid->InsertNextCell( tri->GetCellType(), tri->GetPointIds() );
    }
985
    
Allen Sanderson's avatar
Allen Sanderson committed
986
987
988
989
990
    tri->Delete();
  }
  else //if( element_dimension == 3 )
  {
    vtkWedge *wedge = vtkWedge::New();
991

Allen Sanderson's avatar
Allen Sanderson committed
992
993
//    mins[0] = maxs[0] = 187219;

994
    for( int i=mins[0]; i<=maxs[0]; i+=strides[0] )
Allen Sanderson's avatar
Allen Sanderson committed
995
    {
996
997
998
999
1000
      unsigned int index = i * nvertices;
      wedge->GetPointIds()->SetId( 0, index   );
      wedge->GetPointIds()->SetId( 1, index+1 );
      wedge->GetPointIds()->SetId( 2, index+2 );
      wedge->GetPointIds()->SetId( 3, index+3 );
For faster browsing, not all history is shown. View entire blame