Method InterpolateCellToNodes

Filters/General/vtkFiniteElementFieldDistributor.cxx:1736–1802

----------------------------------------------------------------------------

//----------------------------------------------------------------------------
void vtkFiniteElementFieldDistributor::vtkInternals::InterpolateCellToNodes(vtkIdType cellId,
  vtkCellArray* newCells, vtkPoints* newPoints, vtkCellData* oldCd, vtkPointData* hCurlFields,
  vtkPointData* hDivFields)
{
  // we will interpolate onto the points found at new point ids. (from cell explosion)
  const vtkIdType* newPts = nullptr;
  vtkIdType newNpts = 0;
  newCells->GetCellAtId(cellId, newNpts, newPts);
  if (this->Vblps.RequiresInitialization(this->RefElement, newNpts))
  {
    auto pCoords = this->GetLagrangePCoords(this->RefElement, newNpts);
    // for all others, need to shift center of element to (0,0,0)
    if (this->RefElement != VTK_TRIANGLE && this->RefElement != VTK_TETRA)
    {
      std::transform(pCoords.cbegin(), pCoords.cend(), pCoords.begin(),
        [](const double& val) -> double { return 2 * (val - 0.5); });
    }
    this->Vblps.Initialize(this->RefElement, pCoords.data(), newNpts);
    this->CachedParametricCoordinates = pCoords;
  }
  this->Vblps.ComputeJacobian(
    this->RefElement, this->CachedParametricCoordinates.data(), newNpts, newPts, newPoints);

  for (const auto& fieldName : this->hCurlSpec().Fields)
  {
    std::vector<double> coeffs = ::GetEdgeAttributes(fieldName, oldCd, cellId);
    if (coeffs.empty())
    {
      continue;
    }
    vtkDataArray* outArr = hCurlFields->GetArray(fieldName.c_str());
    const auto refVblpMat = this->Vblps.GetVblp(::SpaceType::HCurl, this->RefElement);
    const auto jacMat = this->Vblps.GetJacobianMatrices(this->RefElement);
    if (refVblpMat != nullptr && jacMat != nullptr)
    {
      const auto physVblpMat =
        this->Vblps.ReferenceToPhysical(*refVblpMat, *jacMat, SpaceType::HCurl);
      ::InterpolateToNodes(physVblpMat, coeffs, newNpts, newPts, outArr);
    }
  }

  for (const auto& fieldName : this->hDivSpec().Fields)
  {
    std::vector<double> coeffs;
    if (this->RefElement == VTK_QUAD || this->RefElement == VTK_TRIANGLE)
    {
      coeffs = ::GetEdgeAttributes(fieldName, oldCd, cellId);
    }
    else
    {
      coeffs = ::GetFaceAttributes(fieldName, oldCd, cellId);
    }
    if (coeffs.empty())
    {
      continue;
    }
    vtkDataArray* outArr = hDivFields->GetArray(fieldName.c_str());
    const auto refVblpMat = this->Vblps.GetVblp(::SpaceType::HDiv, this->RefElement);