vtkDGInvokeOperator.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
// SPDX-License-Identifier: BSD-3-Clause
#ifndef vtkDGInvokeOperator_h
#define vtkDGInvokeOperator_h
 
#include "vtkCellAttribute.h"         // For CellTypeInfo.
#include "vtkDGOperatorEntry.h"       // For ivars.
#include "vtkDataArray.h"             // For implementation.
#include "vtkFiltersCellGridModule.h" // For export macro.
 
#include <array>
#include <cstring> // For std::memset.
#include <vector>
 
VTK_ABI_NAMESPACE_BEGIN
 
class VTKFILTERSCELLGRID_EXPORT vtkDGInvokeOperator
{
public:
  vtkDGInvokeOperator() = default;
  vtkDGInvokeOperator(const vtkDGInvokeOperator&) = default;
 
  struct FetchUnsharedCellDOF
  {
    vtkDataArray* Coefficients{ nullptr };
 
    FetchUnsharedCellDOF() = default;
    FetchUnsharedCellDOF(vtkDataArray* vals, vtkDataArray* conn) { this->Initialize(vals, conn); }
 
    void Initialize(vtkDataArray* vals, vtkDataArray* conn = nullptr)
    {
      (void)conn;
      this->Coefficients = vals;
    }
 
    void operator()(vtkIdType ii, std::vector<double>& tuple)
    {
      // tuple.resize(this->Coefficients->GetNumberOfComponents());
      this->Coefficients->GetTuple(ii, tuple.data());
    }
  };
 
  struct FetchSharedCellDOF
  {
    vtkDataArray* Coefficients{ nullptr };
    vtkDataArray* Connectivity{ nullptr };
    vtkIdType Stride{ 0 };
    std::vector<vtkTypeInt64> ConnTuple;
 
    FetchSharedCellDOF() = default;
    FetchSharedCellDOF(vtkDataArray* vals, vtkDataArray* conn) { this->Initialize(vals, conn); }
 
    void Initialize(vtkDataArray* vals, vtkDataArray* conn)
    {
      this->Coefficients = vals;
      this->Connectivity = conn;
      this->Stride = vals->GetNumberOfComponents();
      this->ConnTuple.resize(conn->GetNumberOfComponents());
    }
 
    void operator()(vtkIdType ii, std::vector<double>& tuple)
    {
      this->Connectivity->GetIntegerTuple(ii, this->ConnTuple.data());
      auto valIt = tuple.begin();
      for (const auto& valId : this->ConnTuple)
      {
        this->Coefficients->GetTuple(valId, &(*valIt));
        valIt += this->Stride;
      }
    }
  };
 
  template <typename InputIterator, typename OutputIterator>
  bool InvokeOp(const vtkDGOperatorEntry& op, const vtkCellAttribute::CellTypeInfo& info,
    InputIterator begin, InputIterator end, OutputIterator out);
 
 
  template <typename InputIterator, typename OutputIterator>
  bool InvokeSharedDOF(const vtkDGOperatorEntry& op, const vtkCellAttribute::CellTypeInfo& info,
    InputIterator begin, InputIterator end, OutputIterator out);
  template <typename InputIterator, typename OutputIterator>
  bool InvokeUnsharedDOF(const vtkDGOperatorEntry& op, const vtkCellAttribute::CellTypeInfo& info,
    InputIterator begin, InputIterator end, OutputIterator out);
 
  template <typename OutputIterator>
  void InnerProduct(int coefficientsPerDOF, int stride, const vtkDGOperatorEntry& op,
    const std::array<double, 3>& rst, OutputIterator& out);
 
  bool Invoke(const vtkDGOperatorEntry& op, const vtkCellAttribute::CellTypeInfo& info,
    std::size_t num, const vtkIdType* cellIds, const double* rst, double* result);
 
  std::vector<double> CoeffTuple;
 
  std::vector<double> OperatorTuple;
 
  FetchSharedCellDOF SharedFetcher;
 
  FetchUnsharedCellDOF DiscontinuousFetcher;
};
 
template <typename InputIterator, typename OutputIterator>
bool vtkDGInvokeOperator::InvokeOp(const vtkDGOperatorEntry& op,
  const vtkCellAttribute::CellTypeInfo& info, InputIterator begin, InputIterator end,
  OutputIterator out)
{
  if (info.DOFSharing.IsValid())
  {
    return this->InvokeSharedDOF(op, info, begin, end, out);
  }
  else
  {
    return this->InvokeUnsharedDOF(op, info, begin, end, out);
  }
}
 
template <typename InputIterator, typename OutputIterator>
bool vtkDGInvokeOperator::InvokeSharedDOF(const vtkDGOperatorEntry& op,
  const vtkCellAttribute::CellTypeInfo& info, InputIterator begin, InputIterator end,
  OutputIterator out)
{
  using namespace vtk::literals;
 
  if (!op)
  {
    return false;
  }
  auto valsit = info.ArraysByRole.find("values"_token);
  if (valsit == info.ArraysByRole.end() || !valsit->second)
  {
    return false;
  }
  auto connit = info.ArraysByRole.find("connectivity"_token);
  if (connit == info.ArraysByRole.end() || !connit->second)
  {
    return false;
  }
  this->SharedFetcher.Initialize(
    vtkDataArray::SafeDownCast(valsit->second), vtkDataArray::SafeDownCast(connit->second));
  int coefficientsPerDOF = valsit->second->GetNumberOfComponents();
  this->CoeffTuple.resize(connit->second->GetNumberOfComponents() * coefficientsPerDOF);
  this->OperatorTuple.resize(op.OperatorSize * op.NumberOfFunctions);
  // There must be exactly one connectivity entry per operator function.
  if (connit->second->GetNumberOfComponents() != op.NumberOfFunctions)
  {
    return false;
  }
  int stride = op.OperatorSize * coefficientsPerDOF;
  vtkIdType lastCellId = -1; // Start with an invalid cell ID.
  for (auto initer = begin; initer != end; ++initer)
  {
    if (lastCellId != initer.GetCellId())
    {
      lastCellId = initer.GetCellId();
      this->SharedFetcher(lastCellId, this->CoeffTuple);
    }
    this->InnerProduct(coefficientsPerDOF, stride, op, initer.GetParameter(), out);
  }
  return true;
}
 
template <typename InputIterator, typename OutputIterator>
bool vtkDGInvokeOperator::InvokeUnsharedDOF(const vtkDGOperatorEntry& op,
  const vtkCellAttribute::CellTypeInfo& info, InputIterator begin, InputIterator end,
  OutputIterator out)
{
  using namespace vtk::literals;
 
  if (!op)
  {
    return false;
  }
  auto valsit = info.ArraysByRole.find("values"_token);
  if (valsit == info.ArraysByRole.end() || !valsit->second)
  {
    return false;
  }
  this->DiscontinuousFetcher.Initialize(vtkDataArray::SafeDownCast(valsit->second));
  int totalCoefficients = valsit->second->GetNumberOfComponents();
  if (totalCoefficients % op.NumberOfFunctions != 0)
  {
    return false;
  }
  int coefficientsPerDOF = totalCoefficients / op.NumberOfFunctions;
  this->CoeffTuple.resize(totalCoefficients);
  this->OperatorTuple.resize(op.OperatorSize * op.NumberOfFunctions);
  int stride = op.OperatorSize * coefficientsPerDOF;
  vtkIdType lastCellId = -1; // Start with an invalid cell ID.
  for (auto initer = begin; initer != end; ++initer)
  {
    if (lastCellId != initer.GetCellId())
    {
      lastCellId = initer.GetCellId();
      this->DiscontinuousFetcher(lastCellId, this->CoeffTuple);
    }
    this->InnerProduct(coefficientsPerDOF, stride, op, initer.GetParameter(), out);
  }
  return true;
}
 
template <typename OutputIterator>
void vtkDGInvokeOperator::InnerProduct(int coefficientsPerDOF, int stride,
  const vtkDGOperatorEntry& op, const std::array<double, 3>& rst, OutputIterator& out)
{
  // Compute the values for \a op at \a rst.
  op.Op(rst, this->OperatorTuple);
  // Zero the output for this input cell-id+(r,s,t)-tuple.
  std::memset(&(out[0]), 0, sizeof(double) * stride);
  // Now sum the "inner product" of the coeffient tuple and operator tuple into the output.
  for (int ii = 0; ii < coefficientsPerDOF; ++ii)
  {
    for (int jj = 0; jj < op.NumberOfFunctions; ++jj)
    {
      for (int kk = 0; kk < op.OperatorSize; ++kk)
      {
        out[ii * op.OperatorSize + kk] += this->CoeffTuple[jj * coefficientsPerDOF + ii] *
          this->OperatorTuple[jj * op.OperatorSize + kk];
      }
    }
  }
  out += stride;
}
 
VTK_ABI_NAMESPACE_END
#endif // vtkDGInvokeOperator_h