avtIVPNIMRODField.C 16.1 KB
Newer Older
1
2
 /*****************************************************************************
*
brugger's avatar
 
brugger committed
3
* Copyright (c) 2000 - 2011, 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.
*
*****************************************************************************/

// ************************************************************************* //
//                             avtIVPNIMRODField.C                            //
// ************************************************************************* //

#include "avtIVPNIMRODField.h"

#include <DebugStream.h>

47
#include <vtkStructuredGrid.h>
48
49
#include <vtkIntArray.h>
#include <vtkFloatArray.h>
50
#include <vtkIdList.h>
51
52
53
54
55
56
57
58
59
60
61
62

#include <math.h>

// ****************************************************************************
//  Method: avtIVPNIMRODField constructor
//
//  Creationist: Allen Sanderson
//  Creation:   20 November 2009
//
// ****************************************************************************

avtIVPNIMRODField::avtIVPNIMRODField( vtkDataSet* dataset, 
63
64
65
66
                                      avtCellLocator* locator ) : 
    avtIVPVTKField( dataset, locator ),
    grid_fourier_series( 0 ), data_fourier_series( 0 ),
    Drad( 2 ), Dtheta( 2 )
67
{
68
69
  vtkStructuredGrid* grid = vtkStructuredGrid::SafeDownCast( dataset );
  
70
  // Pick off all of the data stored with the vtk field.
71
72
73
  if( grid->GetDataDimension() == 3 )
  {
    int dims[3];
74

75
    grid->GetDimensions( dims );
76

77
78
79
80
    Nrad   = dims[0];
    Ntheta = dims[1];
    Nphi   = dims[2];
  }
81
  else
82
83
84
85
86
  {
    debug1 << "Mesh dimension is not 3 " << endl;
    
    return;
  }
87
88

  // Dummy variable to the template class
89
  float fltVar = 0;
90
91

  // The mesh elements.
92
93
//  grid_fourier_series =
//    SetDataPointer( ds, fltVar, "hidden/grid_fourier_series", 3 );
94
95

  // Vector values from the field.
96
97
98
//  data_fourier_series =
//    SetDataPointer( ds, fltVar, "hidden/data_fourier_series", 3 );

99
100
101
102
103
104
105
106
107
108
}

// ****************************************************************************
//  Method: avtIVPNIMRODField constructor
//
//  Creationist: Allen Sanderson
//  Creation:   20 November 2009
//
// ****************************************************************************

109
110
111
112
113
114
115
116
117
avtIVPNIMRODField::avtIVPNIMRODField( unsigned int nRad,
                                      unsigned int nTheta,
                                      unsigned int nPhi,
                                      float *gfs,
                                      float *dfs ) 
  : avtIVPVTKField(NULL, NULL),
    grid_fourier_series( gfs ), data_fourier_series( dfs ),
    Nrad( nRad ), Ntheta( nTheta ), Nphi( nPhi ),
    Drad( 2 ), Dtheta( 2 )
118
{
119
120
121
122
123
124
125
126
127
128
129
130
131
132
  lagrange_nodes[0][0] = 0.0;
  lagrange_nodes[1][0] = 0.0; lagrange_nodes[1][1] = 1.0;
  lagrange_nodes[2][0] = 0.0; lagrange_nodes[2][1] = 0.5; lagrange_nodes[2][2] = 0.5;
  lagrange_nodes[3][0] = 0.0; lagrange_nodes[3][1] = 0.276393202250021031; lagrange_nodes[3][2] = 0.723606797749978969; lagrange_nodes[3][3] = 1.0;
  lagrange_nodes[4][0] = 0.0; lagrange_nodes[4][1] = 0.172673164646011429; lagrange_nodes[4][2] = 0.500000000000000000; lagrange_nodes[4][3] = 0.827326835353988571; lagrange_nodes[4][4] = 1.0;
  lagrange_nodes[5][0] = 0.0; lagrange_nodes[5][1] = 0.117472338035267654; lagrange_nodes[5][2] = 0.357384241759677453; lagrange_nodes[5][3] = 0.642615758240322547; lagrange_nodes[5][4] = 0.882527661964732346; lagrange_nodes[5][5] = 1.0;

//       { 0.0 },
//       { 0.0, 1.0 },
//       { 0.0, 0.5, 1.0 },
//       { 0.0, 0.276393202250021031, 0.723606797749978969, 1.0 },
//       { 0.0, 0.172673164646011429, 0.500000000000000000, 0.827326835353988571, 1.0 },
//       { 0.0, 0.117472338035267654, 0.357384241759677453, 0.642615758240322547, 0.882527661964732346, 1.0 },
//     };
133
134
135
136
137
138
139
140
141
142
143
144
145
}


// ****************************************************************************
//  Method: avtIVPNIMRODField destructor
//
//  Creationist: Allen Sanderson
//  Creation:   20 November 2009
//
// ****************************************************************************

avtIVPNIMRODField::~avtIVPNIMRODField()
{
146
147
  if( grid_fourier_series ) free(grid_fourier_series);
  if( data_fourier_series ) free(data_fourier_series);
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
}


// ****************************************************************************
//  Method: avtIVPNIMRODField destructor
//
//  Creationist: Allen Sanderson
//  Creation:   20 November 2009
//
// ****************************************************************************

template< class type >
type* avtIVPNIMRODField::SetDataPointer( vtkDataSet *ds,
                                        const type var,
                                        const char* varname,
                                        const int ncomponents )
{
  vtkDataArray *array;
  int XX;

  // Because the triangluar mesh is defined by using non unique points
  // and the data is cell centered data VisIt moves it out to the
  // nodes for STREAMLINES thus there are 3 times the number of
  // original values.
  if( ds->GetPointData()->GetArray(varname) )
  {
    array = ds->GetPointData()->GetArray(varname);
  }

  if( array == 0 )
  {
    debug1 << "Variable " << varname
           << " does not exist"
           << endl;
    
    return 0;
  }

186
187
188
  int ntuples = Nrad*Ntheta*Nphi;

  if( ntuples != array->GetNumberOfTuples() ||
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
      ncomponents != array->GetNumberOfComponents() )
  {
    debug1 << "Variable " << varname
           << " size does not equal the number elements and/or components"
           << endl;
    
    return 0;
  }

  type* newptr = new type[ntuples*ncomponents];

  if( newptr == 0 )
  {
    debug1 << "Variable " << varname << " can not allocated" << endl;
    return 0;
  }

  if( array->IsA("vtkIntArray") ) 
  {
    int* ptr = (int*) array->GetVoidPointer(0);

    for( int i=0; i<ntuples; ++i )
      for( int j=0; j<ncomponents; ++j )
212
        newptr[i*ncomponents+j] = ptr[i*ncomponents+j];
213
214
215
216
217
218
219
220
221

    return newptr;
  }
  else if( array->IsA("vtkFloatArray") ) 
  {
    float* ptr = (float*) array->GetVoidPointer(0);

    for( int i=0; i<ntuples; ++i )
      for( int j=0; j<ncomponents; ++j )
222
        newptr[i*ncomponents+j] = ptr[i*ncomponents+j];
223
224
225
226
227
228
229
230
231

    return newptr;
  }
  else if( array->IsA("vtkDoubleArray") ) 
  {
    double* ptr = (double*) array->GetVoidPointer(0);

    for( int i=0; i<ntuples; ++i )
      for( int j=0; j<ncomponents; ++j )
232
        newptr[i*ncomponents+j] = ptr[i*ncomponents+j];
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

    return newptr;
  }
  else
  {
    debug1 << "avtIVPNIMRODField::SetDataPointer "
           << "Variable " << varname
           << " is not of type float - can not safely down cast"
           << endl;
    
    return 0;
  }
}


// ****************************************************************************
//  Method: avtIVPNIMRODField::operator
//
//  Purpose: Evaluates a point location by consulting a M3D C1 grid.
//      Gets the B field components directly - should not be used for
//      calculating integral curves.
//
//  Programmer: Allen Sanderson
//  Creation:   October 24, 2009
//
//  Modifications:
//
// ****************************************************************************

avtVector
avtIVPNIMRODField::operator()( const double &t, const avtVector &p ) const
{
265
  avtVector linear_vec = avtIVPVTKField::operator()( t, p );
266

267
268
  if( !FindCell( t, p ) )
    throw Undefined();
269

270
  vtkIdList *ptIds = vtkIdList::New();
271

272
  ds->GetCellPoints( lastCell, ptIds );
273

274
  double center[3], pt[3];
275

276
277
278
279
280
281
  for( unsigned int i=0; i<Nphi; ++ i )
  { 
    ds->GetPoint( i+62*101, center);
    cerr << i << "  "
         << center[0] << "  " << center[1] << "  " << center[2] << "   "
         << sqrt(center[0]*center[0]+center[1]*center[1]) << "  " << endl;
282
283
  }

284
  cerr << lastCell << " ( " << endl;
285

286
  unsigned int iRad=Nrad, iTheta=Ntheta, iPhi=Nphi;
287

288
289
  double minRad=+1e10, minTheta=+1e10, minPhi=+1e10;
  double maxRad=-1e10, maxTheta=-1e10, maxPhi=-1e10;
290

291
  unsigned int ic=0;
292
293


294
295
296
  for( unsigned int i=0; i<ptIds->GetNumberOfIds(); ++i )
  {
    unsigned int tRad, tTheta, tPhi, tmp = ptIds->GetId(i);
297

298
    cerr << tmp << " ";
299

300
    tRad = tmp / (Ntheta*Nphi);
301

302
    tmp = tmp % (Ntheta*Nphi);
303

304
    tTheta = tmp / Nphi;
305

306
307
308
309
310
    tPhi = tmp % Nrad;
        
    if( iRad >= tRad )
    {
      iRad = tRad;
311

312
313
314
      if( iTheta >= tTheta )
      {
        iTheta = tTheta;
315

316
317
318
319
320
321
322
323
        if( iPhi > tPhi )
        {
          iPhi = tPhi;
          
          ic = i;
        }
      }
    }
324
325


326
327
    ds->GetPoint( ptIds->GetId(i), pt);
    ds->GetPoint( tPhi, center);
328

329
330
331
      cerr << "  " << tPhi << "  " << tTheta << "  " << tRad << "  "
//       << center[0] << "  " << center[1] << "  " << center[2] << "   ";
         << pt[0] << "  " << pt[1] << "  " << pt[2] << "   ";
332

333
334
    double dx = pt[0] - center[0];
    double dz = pt[2] - center[2];
335

336
337
338
    double rad  = sqrt( dx*dx + dz*dz );
    double theta = atan2( dz, dx );
    double phi   = atan2( pt[1], pt[0] );
339

340
    cerr << phi << " " << theta << " " << rad << " " << endl;
341

342
343
    if( minRad > rad ) minRad = rad;
    if( maxRad < rad ) maxRad = rad;
344

345
346
    if( minTheta > theta ) minTheta = theta;
    if( maxTheta < theta ) maxTheta = theta;
347

348
349
    if( minPhi > phi ) minPhi = phi;
    if( maxPhi < phi ) maxPhi = phi;
350
351
352

  }

353
  cerr << ") " << iPhi << "  " << iTheta << "  " << iRad << endl;
354

355
356
  cerr << minRad << " " << minTheta << " " << minPhi << " " << endl;
  cerr << maxRad << " " << maxTheta << " " << maxPhi << " " << endl;
357
358


359
360
  double x[3];
  ds->GetPoint( ptIds->GetId(ic), x);
361

362
  double rad, theta, phi;
363
364

  
365
366
  vec3 P;
  mat3 DRV;
367

368
//  interpolate( rad, theta, phi, &P, &DRV );
369

370
  return linear_vec;
371
372
373
}

// ****************************************************************************
374
//  Method: avtIVPNIMRODField::ConvertToCartesian
375
//
376
377
//  Purpose:
//      Converts the coordinates from cylindrical to cartesian coordinates
378
//
379
380
//  Programmer: Christoph Garth
//  Creation:   February 25, 2008
381
382
383
//
// ****************************************************************************

384
385
386
387
avtVector 
avtIVPNIMRODField::ConvertToCartesian(const avtVector& pt) const
{
  return pt;
388
389
390
}

// ****************************************************************************
391
//  Method: avtIVPNIMRODField::ConvertToCylindrical
392
//
393
394
395
396
397
//  Purpose:
//      Converts the coordinates from cylindrical to cartesian coordinates
//
//  Programmer: Christoph Garth
//  Creation:   February 25, 2008
398
399
400
//
// ****************************************************************************

401
402
403
404
avtVector 
avtIVPNIMRODField::ConvertToCylindrical(const avtVector& pt) const
{
  return pt;
405
406
}

407
// -------------------------------------------------------------------------
408
409


410
411
412
413
414
415
void avtIVPNIMRODField::lagrange_weights( unsigned int DEG, const double s, 
                                          double* w, double *d ) const
{
    for( int i=0; i<=DEG; ++i )
    {
        double nom = 1.0, den = 1.0, dnom = 0.0;
416

417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
        for( int j=0; j<=DEG; ++j )
        {
            if( i==j )
                continue;

            den *= (lagrange_nodes[DEG][i]-lagrange_nodes[DEG][j]);
            nom *= (s-lagrange_nodes[DEG][j]);
            
            double dtmp = 1.0;
                    
            for( int k=0; k<=DEG; ++k )
            {
                if( k==i || k==j )
                    continue;
            
                dtmp *= (s-lagrange_nodes[DEG][k]);
            }
            
            dnom += dtmp;
        }
437

438
439
440
        w[i] = nom/den;
        d[i] = dnom/den;
    }
441
442
}

443
// -------------------------------------------------------------------------
444

445
446
void avtIVPNIMRODField::fourier_weights( unsigned int N, const double t, 
                                         double* w, double* d ) const
447
{
448
449
450
451
452
453
454
455
456
457
458
459
460
461
    // 0th coefficient
    w[0] = 1.0;
    d[0] = 0.0;
     
    for( int n=1, m=N-1; n<N/2; n++, m-- ) 
    {
        double alpha = 2*M_PI*n;
        
        w[n] =  2.0*cos( alpha*t );
        d[n] = -2.0*sin( alpha*t ) * alpha;
        
        w[m] = -2.0*sin( alpha*t );
        d[m] = -2.0*cos( alpha*t ) * alpha;
    }
462

463
464
465
    // N/2 coefficient
    w[N/2] =  cos( M_PI*t );
    d[N/2] = -sin( M_PI*t ) * M_PI;
466
467
}

468
// -------------------------------------------------------------------------
469

470
471
void avtIVPNIMRODField::interpolate( double rad, double theta, double phi,
                                     vec3* P, mat3* DRV ) const
472
{
473
    float *vecs = data_fourier_series;
474

475
    // rad, theta, phi come in parametrized on a unit space cube.
476

477
    // P is the point in physical space.
478

479
480
481
    // Transform from the unit space cube 0->1 into cell index space.
    rad   *= (Nrad-1)/Drad;
    theta *= (Ntheta-1)/Dtheta;
482

483
484
485
486
487
    // Get the integer offset values for the cell index space. 
    unsigned int qrad =
      std::max( 0.0, std::min( floor(rad), (double)(Nrad-1)/Drad - 1 ) );
    unsigned int qtheta =
      std::max( 0.0, std::min( floor(theta), (double)(Ntheta-1)/Dtheta - 1 ) );
488

489
490
491
    // Subtract the integer offset to get the relative cell index values.
    rad -= qrad;
    theta -= qtheta;
492

493
494
495
    // Index space into original grid.
    qrad *= Drad;
    qtheta *= Dtheta;
496

497
    // ---
498

499
500
    double wtheta[Dtheta+1], dtheta[Dtheta+1];
    lagrange_weights( Dtheta, theta, wtheta, dtheta );
501

502
503
    double wrad[Drad+1], drad[Drad+1];
    lagrange_weights( Drad, rad, wrad, drad );
504

505
506
    double wphi[Nphi-1], dphi[Nphi-1];
    fourier_weights( Nphi-1, phi, wphi, dphi );
507

508
509
510
    // Returned vector value in physical space.
    vec3 p;     // = vec3::Zero();
    vec3 Dp[3]; // = mat3::Zero();
511

512
513
    // vec3 v  = vec3::Zero();
    // mat3 Dv = mat3::Zero();
514

515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
    for( int i=0; i<=Drad; ++i )
    {
        for( int j=0; j<=Dtheta; ++j )
        {
            for( int k=0; k<=Nphi-1; ++k )
            {
              vec3 data =
                vec3( vecs +
                      3 * (i + qrad + (qtheta+j)*Nrad + k*Nrad*Ntheta) );
                                
                p += wphi[k] * wtheta[j] * wrad[i] * data;
                Dp[0] += wphi[k] * wtheta[j] * drad[i] * data * (Nrad-1)/Drad;
                Dp[1] += wphi[k] * dtheta[j] * wrad[i] * data * (Ntheta-1)/Dtheta;
                Dp[2] += dphi[k] * wtheta[j] * wrad[i] * data;
            }

            // data = vecs + qrad + (qtheta+j)*Nrad + k*Nrad*Ntheta;
            // 
            // for( int i=0; i<=Drad; ++i )
            // {
            //     v += wphi[k] * wtheta[j] * wrad[i] * data[i];
            //     Dv.row(0) += wphi[k] * wtheta[j] * drad[i] * data[i] * (Nrad-1)/Drad;
            //     Dv.row(1) += wphi[k] * dtheta[j] * wrad[i] * data[i] * (Ntheta-1)/Dtheta;
            //     Dv.row(2) += dphi[k] * wtheta[j] * drad[i] * data[i];
            // }
        }
    }
    
    *P   = p;
//    *DRV = Dp;
545
}