vtkVisItInterpolatedVelocityField.C 8.06 KB
Newer Older
1 2
/*****************************************************************************
*
3
* Copyright (c) 2000 - 2014, Lawrence Livermore National Security, LLC
4
* Produced at the Lawrence Livermore National Laboratory
5
* LLNL-CODE-442911
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 47 48 49
* 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.
*
*****************************************************************************/

#include <vtkVisItInterpolatedVelocityField.h>

#include <float.h>

#include <vtkCell.h>
#include <vtkDataArray.h>
#include <vtkDataSet.h>
#include <vtkGenericCell.h>
#include <vtkIdList.h>
#include <vtkObjectFactory.h>
#include <vtkPointData.h>
50
#include <vtkCellData.h>
51
#include <vtkRectilinearGrid.h>
52
#include <vtkVisItCellLocator.h>
53 54 55 56

#include <vtkVisItUtility.h>
#include <DebugStream.h>

57 58
static void
InterpVector(vtkGenericCell *cell, int numPts, vtkDataArray *vectors, double *weights, double *vel);
59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79

vtkVisItInterpolatedVelocityField* vtkVisItInterpolatedVelocityField::New()
{
  // First try to create the object from the vtkObjectFactory
  vtkObject* ret = vtkObjectFactory::CreateInstance("vtkVisItInterpolatedVelocityField");
  if(ret)
    {
    return (vtkVisItInterpolatedVelocityField*)ret;
    }
  // If the factory was unable to create the object, then create it here.
  return new vtkVisItInterpolatedVelocityField;
}

vtkVisItInterpolatedVelocityField::vtkVisItInterpolatedVelocityField()
{
   ds = NULL;
   lastCell = 0;
   pcoords[0] = pcoords[1] = pcoords[2] = 0.;
   for (int i = 0 ; i < 1024 ; i++)
      weights[i] = 0.;
   locator = NULL;
80 81 82 83
   doPathlines = false;
   nextTimeName = "";
   curTime = 0.;
   nextTime = 1.;
84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107
}

vtkVisItInterpolatedVelocityField::~vtkVisItInterpolatedVelocityField()
{
   if (ds != NULL)
       ds->Delete();
   if (locator != NULL)
       locator->Delete();
}

void
vtkVisItInterpolatedVelocityField::SetDataSet(vtkDataSet *ds_)
{
   if (ds != NULL)
       ds->Delete();
   if (locator != NULL)
   {
       locator->Delete();
       locator = NULL;
   }
   ds = ds_;
   ds->Register(NULL);
}

108 109 110 111 112 113 114 115 116 117 118 119 120
void
vtkVisItInterpolatedVelocityField::SetLocator(vtkVisItCellLocator *loc)
{
   if (locator != NULL)
   {
       locator->Delete();
       locator = NULL;
   }
   locator = loc;
   if (locator != NULL)
       locator->Register(NULL);
}

121
bool
122
vtkVisItInterpolatedVelocityField::Evaluate(double *pt, double *vel, double t)
123 124 125 126 127 128
{
    if (ds == NULL)
    {
        debug1 <<" vtkVisItInterpolatedVelocityField::No data set to evaluate!" << endl;
        return false;
    }
129 130

    bool nodeCenteredVector = true;
131
    vtkDataArray *vectors =  ds->GetPointData()->GetVectors();
132 133
    if (vectors == NULL)
    {
134 135 136 137 138 139 140
        vectors = ds->GetCellData()->GetVectors();
        if (vectors == NULL)
        {
            debug1 <<" vtkVisItInterpolatedVelocityField::Can't locate vectors to interpolate" << endl;
            return false;
        }
        nodeCenteredVector = false;
141
    }
142 143 144

    vtkDataArray *vectors2 =  NULL;
    if (doPathlines)
145
    {
146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163
        if (nodeCenteredVector)
        {
            vectors2 = ds->GetPointData()->GetArray(nextTimeName.c_str());
        }
        else
        {
            vectors2 = ds->GetCellData()->GetArray(nextTimeName.c_str());
        }

        if (vectors2 == NULL)
        {
            debug1 <<" vtkVisItInterpolatedVelocityField::Can't locate vectors to interpolate" << endl;
            return false;
        }
        if (vectors == vectors2)
        {
            debug1 << "vtkVisItInterpolatedVelocityField::Evaluate - Problem: The two vector fields are the same." << endl;
        }
164 165
    }
    
166
    vtkIdType cell = -1;
167 168 169 170 171 172 173 174

    // This is vtkVisItUtility::FindCell, except we cache the locator.
    // I should probably refactor that method, but I'm short on time.
    if (ds->GetDataObjectType() == VTK_RECTILINEAR_GRID)
    {
        int ijk[3];
        vtkRectilinearGrid *rgrid = (vtkRectilinearGrid*)ds;
        if (vtkVisItUtility::ComputeStructuredCoordinates(rgrid, pt, ijk) == 0)
175
            return false;
176 177 178 179 180 181 182 183
        cell = rgrid->ComputeCellId(ijk);
        if (cell < 0)
            return false;
    }
    else
    {
        if (locator == NULL)
        {
184
            locator = vtkVisItCellLocator::New();
185
            locator->SetDataSet(ds);
186
            locator->IgnoreGhostsOn();
187 188 189
            locator->BuildLocator();
        }

190 191 192
        double rad = 1e-6, dist=0.0;
        double resPt[3]={0.0,0.0,0.0};
        int subId = 0;
pugmire's avatar
pugmire committed
193
        locator->IgnoreGhostsOff();
194 195
        locator->FindClosestPointWithinRadius(pt, rad, resPt,
                                              cell, subId, dist);
196 197 198 199 200 201 202
    }
   
    if (cell < 0)
        return false;

    lastCell = cell;

203 204
    //For zone centered vector fields:
    if (!nodeCenteredVector)
205
    {
206 207 208 209 210 211 212 213 214 215 216
        vectors->GetTuple(cell, vel);
        
        if (doPathlines)
        {
            double vel2[3];
            vectors2->GetTuple(cell, vel2);
            double prop1 = 1. - (t - curTime) / (nextTime - curTime);
            vel[0] = prop1*vel[0] + (1-prop1)*vel2[0];
            vel[1] = prop1*vel[1] + (1-prop1)*vel2[1];
            vel[2] = prop1*vel[2] + (1-prop1)*vel2[2];
        }
217
    }
218
    else
219
    {
220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237
        vtkGenericCell *GenCell = vtkGenericCell::New();
        ds->GetCell(cell, GenCell);
        
        int numPts = GenCell->GetNumberOfPoints();

        double closestPoint[3], dist2;
        int subId;
        int val = GenCell->EvaluatePosition(pt, closestPoint, subId, pcoords, dist2, weights);

        if (val <= 0)
        {
            GenCell->Delete();
            return false;
        }
        // interpolate the vectors
        InterpVector(GenCell, numPts, vectors, weights, vel);

        if (doPathlines)
238
        {
239 240 241 242 243 244 245
            double vel2[3];
            InterpVector(GenCell, numPts, vectors2, weights, vel2);
            
            double prop1 = 1. - (t - curTime) / (nextTime - curTime);
            vel[0] = prop1*vel[0] + (1-prop1)*vel2[0];
            vel[1] = prop1*vel[1] + (1-prop1)*vel2[1];
            vel[2] = prop1*vel[2] + (1-prop1)*vel2[2];
246
        }
247
        GenCell->Delete();
248 249 250 251
    }

    return true;
}
252 253 254 255 256 257 258 259 260 261 262 263 264 265 266


static void
InterpVector(vtkGenericCell *cell, int numPts, vtkDataArray *vectors, double *weights, double *vel)
{
    vel[0] = vel[1] = vel[2] = 0;
    double vec[3];
    for (int j=0; j < numPts; j++)
    {
        int id = cell->PointIds->GetId(j);
        vectors->GetTuple(id, vec);
        for (int i=0; i < 3; i++)
            vel[i] += vec[i] * weights[j];
    }
}