------------------------------------------------------------------------------
| 1265 | |
| 1266 | //------------------------------------------------------------------------------ |
| 1267 | double vtkLagrangianParticleTracker::ComputeCellLength(vtkLagrangianParticle* particle) |
| 1268 | { |
| 1269 | double cellLength = 1.0; |
| 1270 | vtkDataSet* dataset = nullptr; |
| 1271 | vtkGenericCell* cell = particle->GetThreadedData()->GenericCell; |
| 1272 | if (!cell) |
| 1273 | { |
| 1274 | vtkErrorMacro("Could not recover a generic cell for cell length computation"); |
| 1275 | return 1.0; |
| 1276 | } |
| 1277 | |
| 1278 | vtkIdType cellId; |
| 1279 | vtkAbstractCellLocator* loc; |
| 1280 | double* weights; |
| 1281 | if (this->IntegrationModel->FindInLocators( |
| 1282 | particle->GetPosition(), particle, dataset, cellId, loc, weights)) |
| 1283 | { |
| 1284 | dataset->GetCell(cellId, cell); |
| 1285 | } |
| 1286 | else |
| 1287 | { |
| 1288 | return -1.0; // no cell found |
| 1289 | } |
| 1290 | |
| 1291 | double* vel = particle->GetVelocity(); |
| 1292 | if (this->CellLengthComputationMode == STEP_CUR_CELL_VEL_DIR && vtkMath::Norm(vel) > 0.0) |
| 1293 | { |
| 1294 | double velHat[3] = { vel[0], vel[1], vel[2] }; |
| 1295 | vtkMath::Normalize(velHat); |
| 1296 | double tmpCellLength = 0.0; |
| 1297 | for (int ne = 0; ne < cell->GetNumberOfEdges(); ++ne) |
| 1298 | { |
| 1299 | double evect[3], x0[3], x1[3]; |
| 1300 | vtkCell* edge = cell->GetEdge(ne); |
| 1301 | vtkIdType e0 = edge->GetPointId(0); |
| 1302 | vtkIdType e1 = edge->GetPointId(1); |
| 1303 | |
| 1304 | dataset->GetPoint(e0, x0); |
| 1305 | dataset->GetPoint(e1, x1); |
| 1306 | vtkMath::Subtract(x0, x1, evect); |
| 1307 | double elength = std::fabs(vtkMath::Dot(evect, velHat)); |
| 1308 | tmpCellLength = std::max(tmpCellLength, elength); |
| 1309 | } |
| 1310 | cellLength = tmpCellLength; |
| 1311 | } |
| 1312 | else if (this->CellLengthComputationMode == STEP_CUR_CELL_DIV_THEO && vtkMath::Norm(vel) > 0.0) |
| 1313 | { |
| 1314 | double velHat[3] = { vel[0], vel[1], vel[2] }; |
| 1315 | vtkMath::Normalize(velHat); |
| 1316 | double xa = 0.0; // cell cross-sectional area in velHat direction |
| 1317 | double vol = 0.0; // cell volume |
| 1318 | for (int nf = 0; nf < cell->GetNumberOfFaces(); ++nf) |
| 1319 | { |
| 1320 | double norm[3]; // cell face normal |
| 1321 | double centroid[3] = { 0.0, 0.0, 0.0 }; // cell face centroid |
| 1322 | vtkCell* face = cell->GetFace(nf); |
| 1323 | vtkPoints* pts = face->GetPoints(); |
| 1324 | vtkIdType nPoints = pts->GetNumberOfPoints(); |
no test coverage detected