GaussCubature Class Reference

Detailed Description

Definition at line 57 of file Cubature.H.

Public Member Functions
	GaussCubature (vtkDataSet *_dataSet)
	GaussCubature (vtkDataSet *_dataSet, const std::vector< int > &arrayIDs)
	~GaussCubature ()
	GaussCubature (const GaussCubature &that)=delete
GaussCubature &	operator= (const GaussCubature &that)=delete
void	constructGaussMesh ()
pntDataPairVec	getGaussPointsAndDataAtCell (int cellID)
	Returns coordinates of gauss points and associated data at cell. More...
void	interpolateToGaussPoints ()
	Interpolate values to gauss points for arrays specified in arrayIDs at construction time. More...
void	interpolateToGaussPoints (const std::vector< std::string > &newArrayNames)
	Interpolate values to gauss points for arrays specified in arrayIDs at construction time and specify new names (in order). More...
std::vector< std::vector< double > >	integrateOverAllCells ()
	Integrate arrays specified by arrayIDs over all cells. More...
std::vector< std::vector< double > >	integrateOverAllCells (const std::vector< std::string > &newArrayNames, bool computeRMSE)
	Same as integrateOverAllCells() but with integral value array names specified by newArrayNames and an optional flag to compute the root mean square error (RMSE) integral over an array with error values, an error metric normalized by cell volume. More...
void	setArrayIDs (const std::vector< int > &_arrayIDs)
double	computeJacobian (vtkSmartPointer< vtkGenericCell > genCell, int cellType) const
double	computeCellVolume (vtkSmartPointer< vtkGenericCell > genCell, int cellType) const
vtkDataSet *	getDataSet () const
vtkSmartPointer< vtkPolyData >	getGaussMesh () const
vtkQuadratureSchemeDefinition **	getDict () const
std::vector< int >	getNumComponents () const
std::vector< int >	getArrayIDs () const
int	getTotalComponents () const
void	writeGaussMesh (const char *name) const

Static Public Member Functions
static GaussCubature *	Create (vtkDataSet *_dataSet)
static GaussCubature *	Create (vtkDataSet *_dataSet, const std::vector< int > &arrayIDs)
static std::unique_ptr< GaussCubature >	CreateUnique (vtkDataSet *_dataSet)
static std::unique_ptr< GaussCubature >	CreateUnique (vtkDataSet *_dataSet, const std::vector< int > &arrayIDs)
static std::shared_ptr< GaussCubature >	CreateShared (vtkDataSet *_dataSet)
static std::shared_ptr< GaussCubature >	CreateShared (vtkDataSet *_dataSet, const std::vector< int > &arrayIDs)

Private Member Functions
int	getOffset (int cellID) const
int	interpolateToGaussPointsAtCell (int cellID, vtkSmartPointer< vtkGenericCell > genCell, const std::vector< vtkSmartPointer< vtkDataArray >> &das, std::vector< vtkSmartPointer< vtkDoubleArray >> &daGausses) const
void	integrateOverCell (int cellID, vtkSmartPointer< vtkGenericCell > genCell, vtkSmartPointer< vtkPointData > pd, std::vector< vtkSmartPointer< vtkDoubleArray >> &integralData, std::vector< std::vector< double >> &totalIntegralData) const
void	integrateOverCell (int cellID, vtkSmartPointer< vtkGenericCell > genCell, vtkSmartPointer< vtkPointData > pd, std::vector< vtkSmartPointer< vtkDoubleArray >> &integralData, std::vector< std::vector< double >> &totalIntegralData, const std::vector< std::string > &newArrayNames, bool normalizeByVol) const

Private Attributes
vtkSmartPointer< vtkDataSet >	dataSet
int	numVolCells
vtkSmartPointer< vtkPolyData >	gaussMesh
vtkQuadratureSchemeDefinition **	dict
std::vector< int >	arrayIDs
std::vector< int >	numComponents
int	totalComponents

Constructor & Destructor Documentation

GaussCubature() [1/3]

GaussCubature::GaussCubature ( vtkDataSet *  _dataSet )

explicit

Definition at line 107 of file Cubature.C.

References constructGaussMesh(), and dataSet.

Referenced by Create().

    : dataSet(_dataSet), numVolCells(0), totalComponents(0)
{
  dataSet->GetCellData()->RemoveArray("QuadratureOffSet");
  constructGaussMesh();
}

GaussCubature() [2/3]

GaussCubature::GaussCubature	(	vtkDataSet *	_dataSet,
		const std::vector< int > &	arrayIDs
	)

Definition at line 114 of file Cubature.C.

References constructGaussMesh(), dataSet, and interpolateToGaussPoints().

    : dataSet(_dataSet), numVolCells(0), arrayIDs(_arrayIDs),
      totalComponents(0)
{
  dataSet->GetCellData()->RemoveArray("QuadratureOffSet");
  constructGaussMesh();
  interpolateToGaussPoints();
}

~GaussCubature()

GaussCubature::~GaussCubature ( )

inline

Definition at line 62 of file Cubature.H.

{ delete[] dict; }

GaussCubature() [3/3]

GaussCubature::GaussCubature ( const GaussCubature &  that )

delete

Member Function Documentation

computeCellVolume()

double GaussCubature::computeCellVolume	(	vtkSmartPointer< vtkGenericCell >	genCell,
		int	cellType
	)		const

Definition at line 301 of file Cubature.C.

Referenced by integrateOverCell().

{
  switch (cellType)
  {
    case VTK_TRIANGLE:
    {
      return vtkMeshQuality::TriangleArea(genCell);
    }
    case VTK_QUAD:
    {
      return vtkMeshQuality::QuadArea(genCell);
    }
    case VTK_TETRA:
    {
      return vtkMeshQuality::TetVolume(genCell);
    }
    case VTK_HEXAHEDRON:
    {
      return vtkMeshQuality::HexVolume(genCell);
    }
    default:
    {
      std::cerr << "Error: Cell type: " << cellType << "found "
                << "with no quadrature definition provided" << std::endl;
      exit(1);
    }
  }
}

computeJacobian()

double GaussCubature::computeJacobian	(	vtkSmartPointer< vtkGenericCell >	genCell,
		int	cellType
	)		const

Definition at line 333 of file Cubature.C.

Referenced by integrateOverCell().

{
  switch (cellType)
  {
    case VTK_TRIANGLE:
    {
      return vtkMeshQuality::TriangleArea(genCell);
    }
    case VTK_QUAD:
    {
      return vtkMeshQuality::QuadArea(genCell) / 4.0;
    }
    case VTK_TETRA:
    {
      return vtkMeshQuality::TetVolume(genCell);
    }
    case VTK_HEXAHEDRON:
    {
      return vtkMeshQuality::HexVolume(genCell) / 8.0;
    }
    default:
    {
      std::cerr << "Error: Cell type: " << cellType << "found "
                << "with no quadrature definition provided" << std::endl;
      exit(1);
    }
  }
}

constructGaussMesh()

void GaussCubature::constructGaussMesh

(

)

Definition at line 166 of file Cubature.C.

References arrayIDs, cellType, dataSet, dict, gaussMesh, HEX8, HEX8W, NEM::MSH::New(), numVolCells, TET4, TET4W, TRI3, and TRI3W.

Referenced by GaussCubature().

{
  // check whether arrayIDs exist in mesh
  {
    vtkSmartPointer<vtkPointData> pd = dataSet->GetPointData();
    int numArr = pd->GetNumberOfArrays();
    for (int arrayID : arrayIDs)
    {
      if (arrayID >= numArr)
      {
        std::cerr << "Array ID " << arrayID << " not found in dataset"
                  << std::endl;
        exit(1);
      }
    }
  }

  // Get the dictionary key     
  auto key = vtkQuadratureSchemeDefinition::DICTIONARY();

  // Get the cell types used by the data set
  auto cellTypes = vtkSmartPointer<vtkCellTypes>::New();
  dataSet->GetCellTypes(cellTypes);
  int nCellTypes = cellTypes->GetNumberOfTypes();

  // create offset array and store the dictionary within
  auto offsets = vtkSmartPointer<vtkIdTypeArray>::New();
  std::string basename = "QuadratureOffset";
  offsets->SetName(basename.c_str());
  auto info = vtkSmartPointer<vtkInformation>::New();
  info = offsets->GetInformation();

  for (int typeId = 0; typeId < nCellTypes; ++typeId)
  {
    int cellType = cellTypes->GetCellType(typeId);
    // Initialize quadrature scheme definition for given cell type
    auto def = vtkSmartPointer<vtkQuadratureSchemeDefinition>::New();
    switch (cellType)
    {
      case VTK_TRIANGLE:
        def->Initialize(VTK_TRIANGLE, 3, 3, TRI3, TRI3W);
        break;
      case VTK_TETRA:
        def->Initialize(VTK_TETRA, 4, 4, TET4, TET4W);
        break;
      case VTK_HEXAHEDRON:
        def->Initialize(VTK_HEXAHEDRON, 8, 8, HEX8, HEX8W);
        break;
      default:
        std::cerr << "Error: Cell type: " << cellType << " found "
                  << "with no quadrature definition provided" << std::endl;
        exit(1);
    }
    // the definition must appear in the dictionary associated with
    // the offset array
    key->Set(info, def, cellType);
  }
  // get dictionary size 
  int dictSize = key->Size(info);
  dict = new vtkQuadratureSchemeDefinition *[dictSize];
  key->GetRange(info, dict, 0, 0, dictSize);
  vtkIdType numCells = dataSet->GetNumberOfCells();

  offsets->SetNumberOfTuples(numCells);

#ifdef HAVE_OPENMP
  std::vector<vtkIdType> chunkOffsets(omp_get_max_threads() + 1, 0);
  #pragma omp parallel default(none) shared(numCells, chunkOffsets, offsets)
  {
    vtkIdType subOffset = 0;
    auto genCell = vtkSmartPointer<vtkGenericCell>::New();

    #pragma omp single
    {
      // GetCell method is thread safe when first called from single thread
      dataSet->GetCell(0, genCell);
    }

    #pragma omp for schedule(static)
    for (vtkIdType cellId = 0; cellId < numCells; ++cellId) {
      offsets->SetValue(cellId, subOffset);

      dataSet->GetCell(cellId, genCell);
      int cellType = genCell->GetCellType();
      if (cellType >= VTK_TETRA) numVolCells += 1;

      vtkQuadratureSchemeDefinition *celldef = dict[cellType];
      subOffset += celldef->GetNumberOfQuadraturePoints();
    }

    int tid = omp_get_thread_num();
    chunkOffsets[tid + 1] = subOffset;
    #pragma omp barrier
    #pragma omp single
    {
      for (int i = 1; i < chunkOffsets.size(); ++i) {
        chunkOffsets[i] = chunkOffsets[i] + chunkOffsets[i - 1];
      }
    }

    #pragma omp for schedule(static)
    for (vtkIdType cellId = 0; cellId < numCells; ++cellId) {
      vtkIdType offset = chunkOffsets[tid] + offsets->GetValue(cellId);
      offsets->SetValue(cellId, offset);
    }
  }
#else
  vtkIdType offset = 0;

  for (int cellid = 0; cellid < numCells; ++cellid) {
    offsets->SetValue(cellid, offset);
    int cellType = dataSet->GetCell(cellid)->GetCellType();
    if (cellType >= VTK_TETRA) {
      numVolCells += 1;
    }
    vtkQuadratureSchemeDefinition *celldef = dict[cellType];
    offset += celldef->GetNumberOfQuadraturePoints();
  }
#endif

  dataSet->GetCellData()->AddArray(offsets);

  auto pointGen = vtkSmartPointer<vtkQuadraturePointsGenerator>::New();

  pointGen->SetInputArrayToProcess
      (0, 0, 0,
       vtkDataObject::FIELD_ASSOCIATION_CELLS,
       "QuadratureOffset");
  pointGen->SetInputData(dataSet);
  gaussMesh = vtkSmartPointer<vtkPolyData>::New();
  gaussMesh = vtkPolyData::SafeDownCast(pointGen->GetOutput());
  pointGen->Update();
}

Create() [1/2]

GaussCubature * GaussCubature::Create ( vtkDataSet *  _dataSet )

static

Definition at line 124 of file Cubature.C.

References GaussCubature().

Referenced by CreateShared(), and CreateUnique().

{
  return new GaussCubature(_dataSet);
}

Create() [2/2]

GaussCubature * GaussCubature::Create ( vtkDataSet *  _dataSet, const std::vector< int > &  arrayIDs  )

static

Definition at line 129 of file Cubature.C.

References GaussCubature().

{
  return new GaussCubature(_dataSet, arrayIDs);
}

CreateShared() [1/2]

std::shared_ptr< GaussCubature > GaussCubature::CreateShared ( vtkDataSet *  _dataSet )

static

Definition at line 149 of file Cubature.C.

References Create().

{
  std::shared_ptr<GaussCubature> cuby;
  cuby.reset(GaussCubature::Create(_dataSet));
  return cuby;
}

CreateShared() [2/2]

std::shared_ptr< GaussCubature > GaussCubature::CreateShared ( vtkDataSet *  _dataSet, const std::vector< int > &  arrayIDs  )

static

Definition at line 157 of file Cubature.C.

References Create().

{
  std::shared_ptr<GaussCubature> cuby;
  cuby.reset(GaussCubature::Create(_dataSet, arrayIDs));
  return cuby;
}

CreateUnique() [1/2]

std::unique_ptr< GaussCubature > GaussCubature::CreateUnique ( vtkDataSet *  _dataSet )

static

Definition at line 135 of file Cubature.C.

References Create().

Referenced by meshBase::integrateOverMesh(), and PatchRecovery::PatchRecovery().

{
  return std::unique_ptr<GaussCubature>(GaussCubature::Create(_dataSet));
}

CreateUnique() [2/2]

std::unique_ptr< GaussCubature > GaussCubature::CreateUnique ( vtkDataSet *  _dataSet, const std::vector< int > &  arrayIDs  )

static

Definition at line 141 of file Cubature.C.

References Create().

{
  return std::unique_ptr<GaussCubature>(
      GaussCubature::Create(_dataSet, arrayIDs));
}

getArrayIDs()

std::vector<int> GaussCubature::getArrayIDs ( ) const

inline

Definition at line 118 of file Cubature.H.

{ return arrayIDs; }

getDataSet()

vtkDataSet* GaussCubature::getDataSet ( ) const

inline

Definition at line 113 of file Cubature.H.

{ return dataSet; }

getDict()

vtkQuadratureSchemeDefinition** GaussCubature::getDict ( ) const

inline

Definition at line 116 of file Cubature.H.

{ return dict; }

getGaussMesh()

vtkSmartPointer<vtkPolyData> GaussCubature::getGaussMesh ( ) const

inline

Definition at line 115 of file Cubature.H.

{ return gaussMesh; }

getGaussPointsAndDataAtCell()

pntDataPairVec GaussCubature::getGaussPointsAndDataAtCell

(

int

cellID

)

Parameters

cellId

cell id

Returns

pntDataPairVec (see alias in Cubature.H)

Definition at line 374 of file Cubature.C.

References arrayIDs, data, dataSet, dict, gaussMesh, getOffset(), interpolateToGaussPoints(), interpolateToGaussPointsAtCell(), numComponents, and totalComponents.

{
  if (arrayIDs.empty())
  {
    std::cerr << "no array have been selected for interpolation" << std::endl;
    exit(1);
  }

  int numDataArr = gaussMesh->GetPointData()->GetNumberOfArrays();

  if (numDataArr == 0)
  {
    interpolateToGaussPoints();
  }

  // get number of gauss points in cell from dictionary
  int numGaussPoints = dict[dataSet->GetCell(
      cellID)->GetCellType()]
      ->GetNumberOfQuadraturePoints();
  // get offset from nodeMesh for lookup of gauss points in polyData
  int offset = getOffset(cellID);

  //pntDataPairVec container;
  pntDataPairVec container(numGaussPoints);

  vtkSmartPointer<vtkPointData> pd = gaussMesh->GetPointData();
  for (int i = 0; i < numGaussPoints; ++i)
  {
    double x_tmp[3];
    gaussMesh->GetPoint(offset + i, x_tmp);
    std::vector<double> gaussPnt(x_tmp, x_tmp + 3);
    std::vector<double> data(totalComponents);
    int currcomp = 0;
    for (int j = 0; j < numComponents.size(); ++j)
    {
      double *comps = new double[numComponents[j]];
      pd->GetArray(j)->GetTuple(offset + i, comps);
      for (int k = 0; k < numComponents[j]; ++k)
      {
        data[currcomp] = comps[k];
        ++currcomp;
      }
      delete[] comps;
    }
    container[i] = std::make_pair(gaussPnt, data);
  }
  return container;
}

getNumComponents()

std::vector<int> GaussCubature::getNumComponents ( ) const

inline

Definition at line 117 of file Cubature.H.

{ return numComponents; }

getOffset()

int GaussCubature::getOffset ( int  cellID ) const

private

Definition at line 364 of file Cubature.C.

References dataSet.

Referenced by getGaussPointsAndDataAtCell(), integrateOverCell(), and interpolateToGaussPointsAtCell().

{
  vtkIdType offsets[1];
  vtkIdTypeArray::FastDownCast(
      dataSet->GetCellData()->GetArray("QuadratureOffset"))
      ->GetTypedTuple(cellID, offsets);
  return offsets[0];
}

getTotalComponents()

int GaussCubature::getTotalComponents ( ) const

inline

Definition at line 119 of file Cubature.H.

{ return totalComponents; }

integrateOverAllCells() [1/2]

std::vector< std::vector< double > > GaussCubature::integrateOverAllCells

(

)

Returns

Integral data (first index - array, second index - component)

Definition at line 697 of file Cubature.C.

References arrayIDs, dataSet, gaussMesh, id, integrateOverCell(), interpolateToGaussPoints(), NEM::MSH::New(), and numComponents.

Referenced by meshBase::integrateOverMesh().

{
  int numCells = dataSet->GetNumberOfCells();

  if (gaussMesh->GetPointData()->GetNumberOfArrays() == 0)
  {
    interpolateToGaussPoints();
  }

  vtkSmartPointer<vtkPointData> pd = gaussMesh->GetPointData();
  std::vector<vtkSmartPointer<vtkDoubleArray>> cellIntegralData(arrayIDs.size());
  std::vector<std::vector<double>> totalIntegralData(arrayIDs.size());
  for (int id = 0; id < arrayIDs.size(); ++id)
  {
    std::string arrName(
        dataSet->GetPointData()->GetArrayName(arrayIDs[id]));
    arrName.append("Integral");
    vtkSmartPointer<vtkDoubleArray> integralDatum = vtkSmartPointer<vtkDoubleArray>::New();
    integralDatum->SetName(&arrName[0u]);
    integralDatum->SetNumberOfComponents(numComponents[id]);
    integralDatum->SetNumberOfTuples(dataSet->GetNumberOfCells());
    cellIntegralData[id] = integralDatum;
    totalIntegralData[id].resize(numComponents[id], 0);
  }

#ifdef HAVE_OPENMP
  #pragma omp parallel default(none) shared(numCells, pd, cellIntegralData, totalIntegralData, cerr)
  {
    auto genCell = vtkSmartPointer<vtkGenericCell>::New();

    #pragma omp single
    {
      // GetCell method is thread safe when first called from single thread
      dataSet->GetCell(0, genCell);
    }
    auto partialIntegralData = totalIntegralData;

    #pragma omp for schedule(static)
    for (int i = 0; i < numCells; ++i) {
      integrateOverCell(i, genCell, pd, cellIntegralData, partialIntegralData);
    }

    #pragma omp critical
    {
      // total integral data has same dimensions as cell integral data ...
      // sum partial integrals
      for (int id = 0; id < partialIntegralData.size(); ++id) {
        for (int k = 0; k < partialIntegralData[id].size(); ++k) {
          totalIntegralData[id][k] += partialIntegralData[id][k];
        }
      }
    }
  }
#else
  auto genCell = vtkSmartPointer<vtkGenericCell>::New();
  for (int i = 0; i < numCells; ++i) {
    integrateOverCell(i, genCell, pd, cellIntegralData, totalIntegralData);
  }
#endif

  for (int id = 0; id < arrayIDs.size(); ++id) {
    dataSet->GetCellData()->AddArray(cellIntegralData[id]);
  }
  return totalIntegralData;
}

integrateOverAllCells() [2/2]

std::vector< std::vector< double > > GaussCubature::integrateOverAllCells	(	const std::vector< std::string > &	newArrayNames,
		bool	computeRMSE
	)

Parameters

newArrayNames	integral value array names
computeRMSE	compute RMSE (see above)

Returns

Integral data (first index - array, second index - component)

Definition at line 765 of file Cubature.C.

References dataSet, gaussMesh, id, integrateOverCell(), interpolateToGaussPoints(), and NEM::MSH::New().

{
  interpolateToGaussPoints(newArrayNames);

  vtkSmartPointer<vtkPointData> pd = gaussMesh->GetPointData();
  std::vector<vtkSmartPointer<vtkDoubleArray>> cellIntegralData(newArrayNames.size());
  std::vector<std::vector<double>> totalIntegralData(newArrayNames.size());
  for (int id = 0; id < newArrayNames.size(); ++id)
  {
    std::string name = newArrayNames[id] + "Integral";
    vtkSmartPointer<vtkDoubleArray> integralDatum = vtkSmartPointer<vtkDoubleArray>::New();
    integralDatum->SetName(&name[0u]);
    int numComponent = pd->GetArray(
        &(newArrayNames[id])[0u])->GetNumberOfComponents();
    integralDatum->SetNumberOfComponents(numComponent);
    integralDatum->SetNumberOfTuples(dataSet->GetNumberOfCells());
    cellIntegralData[id] = integralDatum;
    totalIntegralData[id].resize(numComponent, 0);
  }

  int numCells = dataSet->GetNumberOfCells();

#ifdef HAVE_OPENMP
  #pragma omp parallel default(none) \
  shared(numCells, pd, cellIntegralData, totalIntegralData, newArrayNames, computeRMSE)
  {
    auto genCell = vtkSmartPointer<vtkGenericCell>::New();
    #pragma omp single
    {
      // GetCell method is thread safe when first called from single thread
      dataSet->GetCell(0, genCell);
    }
    auto partialIntegralData = totalIntegralData;
    #pragma omp for schedule(static)
    for (int i = 0; i < numCells; ++i) {
      integrateOverCell(i, genCell, pd, cellIntegralData, partialIntegralData,
                        newArrayNames, computeRMSE);
    }
    #pragma omp critical
    {
      // total integral data has same dimensions as cell integral data ...
      // sum partial integrals
      for(int id = 0; id < partialIntegralData.size(); ++id) {
        for(int k = 0; k < partialIntegralData[id].size(); ++k) {
          totalIntegralData[id][k] += partialIntegralData[id][k];
        }
      }
    }
  }
#else
  vtkSmartPointer<vtkGenericCell> genCell = vtkSmartPointer<vtkGenericCell>::New();
  for (int i = 0; i < numCells; ++i) {
    integrateOverCell(i, genCell, pd, cellIntegralData, totalIntegralData,
                      newArrayNames, computeRMSE);
  }
#endif

  for (int id = 0; id < newArrayNames.size(); ++id)
  {
    dataSet->GetCellData()->AddArray(cellIntegralData[id]);
  }
  return totalIntegralData;
}

integrateOverCell() [1/2]

void GaussCubature::integrateOverCell ( int  cellID, vtkSmartPointer< vtkGenericCell >  genCell, vtkSmartPointer< vtkPointData >  pd, std::vector< vtkSmartPointer< vtkDoubleArray >> &  integralData, std::vector< std::vector< double >> &  totalIntegralData  ) const

private

Definition at line 586 of file Cubature.C.

References arrayIDs, cellType, computeJacobian(), data, dataSet, dict, and getOffset().

Referenced by integrateOverAllCells(), and interpolateToGaussPoints().

{
  // putting cell from nodeMesh into genCell
  dataSet->GetCell(cellID, genCell);
  // getting cellType for looking up numGaussPoints in dictionary
  // as well as computing scaled Jacobian
  int cellType = genCell->GetCellType();

  // get number of gauss points in cell from dictionary
  int numGaussPoints;
  const double* quadWeights;
  numGaussPoints = dict[cellType]->GetNumberOfQuadraturePoints();
  quadWeights = dict[cellType]->GetQuadratureWeights();

  // computing Jacobian for integration
  double jacobian = computeJacobian(genCell, cellType);
  // get quadrature weights for this cell type
  // get offset from nodeMesh for lookup of gauss points in polyData
  int offset = getOffset(cellID);
  // holds integrated data for each array
  std::vector<std::vector<double>> data(arrayIDs.size());
  // integration loop
  for (int j = 0; j < integralData.size(); ++j) {
    int numComponent = integralData[j]->GetNumberOfComponents();
    data[j].resize(numComponent, 0.0);
    auto comps = new double[numComponent];
    for (int i = 0; i < numGaussPoints; ++i) {
      pd->GetArray(j)->GetTuple(offset + i, comps);
      for (int k = 0; k < numComponent; ++k) {
        // TODO: generalize to support surface integration
        if (genCell->GetCellDimension() == 3) {
          data[j][k] += comps[k] * quadWeights[i];//*jacobian;
        } else
          data[j][k] += 0.0;
      }
    }
    delete[] comps;
    for (int k = 0; k < numComponent; ++k) {
      data[j][k] *= jacobian;
      totalIntegralData[j][k] += data[j][k];
    }
    // adding integrated value to data of cell
    integralData[j]->SetTuple(cellID, data[j].data());
  }
}

integrateOverCell() [2/2]

void GaussCubature::integrateOverCell ( int  cellID, vtkSmartPointer< vtkGenericCell >  genCell, vtkSmartPointer< vtkPointData >  pd, std::vector< vtkSmartPointer< vtkDoubleArray >> &  integralData, std::vector< std::vector< double >> &  totalIntegralData, const std::vector< std::string > &  newArrayNames, bool  normalizeByVol  ) const

private

Definition at line 638 of file Cubature.C.

References arrayIDs, cellType, computeCellVolume(), computeJacobian(), data, dataSet, dict, and getOffset().

{
  // putting cell from nodeMesh into genCell
  dataSet->GetCell(cellID, genCell);
  // getting cellType for looking up numGaussPoints in dictionary
  int cellType = dataSet->GetCell(cellID)->GetCellType();
  // get number of gauss points in cell from dictionary
  int numGaussPoints = dict[cellType]->GetNumberOfQuadraturePoints();
  // computing Jacobian for integration
  double jacobian = computeJacobian(genCell, cellType);
  // computing volume for RMSE
  double volume = computeCellVolume(genCell, cellType);
  // get quadrature weights for this cell type
  const double *quadWeights = dict[cellType]->GetQuadratureWeights();
  // get offset from nodeMesh for lookup of gauss points in polyData
  int offset = getOffset(cellID);
  // holds integrated data for each array
  std::vector<std::vector<double>> data(arrayIDs.size());
  // integration loop
  for (int j = 0; j < integralData.size(); ++j)
  {
    int numComponent = integralData[j]->GetNumberOfComponents();
    data[j].resize(numComponent);
    auto *comps = new double[numComponent];
    for (int i = 0; i < numGaussPoints; ++i)
    {
      pd->GetArray(&(newArrayNames[j])[0u])->GetTuple(offset + i, comps);
      for (int k = 0; k < numComponent; ++k)
      {
        // TODO: generalize to support surface integration
        if (genCell->GetCellDimension() == 3)
        {
          data[j][k] += comps[k] * quadWeights[i];
        }
        else
          data[j][k] += 0.0;
      }
    }
    delete[] comps;
    // taking sqrt of integrated value (RMSE)
    for (int k = 0; k < numComponent; ++k)
    {
      data[j][k] *= jacobian;
      data[j][k] = (computeRMSE ? std::sqrt(data[j][k] / volume) : data[j][k]);
      //if (computeRMSE)
      //  data[j][k] = std::sqrt(data[j][k] / volume);
      totalIntegralData[j][k] += data[j][k];
    }
    // adding integrated value to data of cell
    integralData[j]->SetTuple(cellID, data[j].data());
  }
}

interpolateToGaussPoints() [1/2]

void GaussCubature::interpolateToGaussPoints

(

)

Definition at line 467 of file Cubature.C.

References arrayIDs, dataSet, gaussMesh, id, interpolateToGaussPointsAtCell(), NEM::MSH::New(), numComponents, and totalComponents.

Referenced by GaussCubature(), getGaussPointsAndDataAtCell(), and integrateOverAllCells().

{
  if (arrayIDs.empty())
  {
    std::cerr << "no arrays selected for interpolation" << std::endl;
    exit(1);
  }

  std::vector<vtkSmartPointer<vtkDoubleArray>> daGausses(arrayIDs.size());
  std::vector<vtkSmartPointer<vtkDataArray>> das(arrayIDs.size());
  numComponents.resize(arrayIDs.size());
  // initializing arrays storing interpolated data
  for (int id = 0; id < arrayIDs.size(); ++id)
  {
    // get desired point data array to be interpolated to gauss points
    vtkSmartPointer<vtkDataArray> da = dataSet->GetPointData()->GetArray(
        arrayIDs[id]);
    // get tuple length of given data
    int numComponent = da->GetNumberOfComponents();
    // declare data array to be populated with values at gauss points
    vtkSmartPointer<vtkDoubleArray> daGauss = vtkSmartPointer<vtkDoubleArray>::New();
    // names and sizing
    daGauss->SetName(
        dataSet->GetPointData()->GetArrayName(arrayIDs[id]));
    daGauss->SetNumberOfComponents(numComponent);
    daGauss->SetNumberOfTuples(gaussMesh->GetNumberOfPoints());
    das[id] = da;
    daGausses[id] = daGauss;
    numComponents[id] = numComponent;
    totalComponents += numComponent;
  }
  int numCells = dataSet->GetNumberOfCells();
#ifdef HAVE_OPENMP
  #pragma omp parallel default(none) shared(numCells, das, daGausses)
  {
    vtkSmartPointer<vtkGenericCell> genCell =
        vtkSmartPointer<vtkGenericCell>::New();
    #pragma omp single
    {
      // GetCell method is thread safe when first called from single thread
      dataSet->GetCell(0, genCell);
    }
    #pragma omp for schedule(static)
    for (int i = 0; i < numCells; ++i) {
      interpolateToGaussPointsAtCell(i, genCell, das, daGausses);
    }
  }
#else
  // generic cell to store given cell in dataSet
  vtkSmartPointer<vtkGenericCell> genCell = vtkSmartPointer<vtkGenericCell>::New();
  for (int i = 0; i < dataSet->GetNumberOfCells(); ++i)
  {
    interpolateToGaussPointsAtCell(i, genCell, das, daGausses);
  }
#endif
  for (int id = 0; id < arrayIDs.size(); ++id)
  {
    gaussMesh->GetPointData()->AddArray(daGausses[id]);
  }
}

interpolateToGaussPoints() [2/2]

void GaussCubature::interpolateToGaussPoints

(

const std::vector< std::string > &

newArrayNames

)

Parameters

newArrayNames

name of transferred data from arrayIDs[i] will be named newArrayNames[i]

Definition at line 529 of file Cubature.C.

References arrayIDs, dataSet, gaussMesh, id, integrateOverCell(), interpolateToGaussPointsAtCell(), and NEM::MSH::New().

                                               {
  if (newArrayNames.empty()) {
    std::cerr << "no arrays selected for interpolation" << std::endl;
    exit(1);
  }

  std::vector<vtkSmartPointer<vtkDoubleArray>> daGausses(newArrayNames.size());
  std::vector<vtkSmartPointer<vtkDataArray>> das(newArrayNames.size());
  // initializing arrays storing interpolated data
  for (int id = 0; id < newArrayNames.size(); ++id) {
    // get desired point data array to be interpolated to gauss points
    vtkSmartPointer<vtkDataArray> da =
        dataSet->GetPointData()->GetArray(&(newArrayNames[id])[0u]);
    // get tuple length of given data
    int numComponent = da->GetNumberOfComponents();
    // declare data array to be populated with values at gauss points
    vtkSmartPointer<vtkDoubleArray> daGauss = vtkSmartPointer<vtkDoubleArray>::New();
    // names and sizing
    daGauss->SetName(&(newArrayNames[id])[0u]);
    daGauss->SetNumberOfComponents(numComponent);
    daGauss->SetNumberOfTuples(gaussMesh->GetNumberOfPoints());
    das[id] = da;
    daGausses[id] = daGauss;
  }
  int numCells = dataSet->GetNumberOfCells();
#ifdef HAVE_OPENMP
  #pragma omp parallel default(none) shared(numCells, das, daGausses)
  {
    // generic cell to store given cell in dataSet
    vtkSmartPointer<vtkGenericCell> genCell =
        vtkSmartPointer<vtkGenericCell>::New();
    #pragma omp single
    {
      // GetCell method is thread safe when first called from single thread
      dataSet->GetCell(0, genCell);
    }
    #pragma omp for schedule(static)
    for (int i = 0; i < numCells; ++i) {
      interpolateToGaussPointsAtCell(i, genCell, das, daGausses);
    }
  }
#else
  // generic cell to store given cell in dataSet
  vtkSmartPointer<vtkGenericCell> genCell =
      vtkSmartPointer<vtkGenericCell>::New();
  for (int i = 0; i < dataSet->GetNumberOfCells(); ++i) {
    interpolateToGaussPointsAtCell(i, genCell, das, daGausses);
  }
#endif
  for (int id = 0; id < arrayIDs.size(); ++id)
  {
    gaussMesh->GetPointData()->AddArray(daGausses[id]);
  }
}

interpolateToGaussPointsAtCell()

int GaussCubature::interpolateToGaussPointsAtCell ( int  cellID, vtkSmartPointer< vtkGenericCell >  genCell, const std::vector< vtkSmartPointer< vtkDataArray >> &  das, std::vector< vtkSmartPointer< vtkDoubleArray >> &  daGausses  ) const

private

Definition at line 425 of file Cubature.C.

References cellType, dataSet, dict, getOffset(), and id.

Referenced by getGaussPointsAndDataAtCell(), and interpolateToGaussPoints().

{
  // putting current cell into genCell
  dataSet->GetCell(cellID, genCell);
  // getting cellType information for lookup in map
  int cellType = dataSet->GetCellType(cellID);
  // get quadrature weights for this cell type
  const double *shapeFunctionWeights = dict[cellType]->GetShapeFunctionWeights();
  // number of gauss points in this cell
  int numGaussPoints = dict[cellType]->GetNumberOfQuadraturePoints();
  // get offset from nodeMesh for lookup of gauss points in polyData
  int offset = getOffset(cellID);
  // interpolation loop
  for (int j = 0; j < numGaussPoints; ++j)
  {
    for (int id = 0; id < das.size(); ++id)
    {
      int numComponent = das[id]->GetNumberOfComponents();
      auto *comps = new double[numComponent];
      std::vector<double> interps(numComponent, 0.0);
      for (int m = 0; m < genCell->GetNumberOfPoints(); ++m)
      {
        int pntId = genCell->GetPointId(m);
        das[id]->GetTuple(pntId, comps);
        for (int h = 0; h < numComponent; ++h)
        {
          interps[h] += comps[h] *
                        shapeFunctionWeights[j * genCell->GetNumberOfPoints() + m];
        }
      }
      delete[] comps;
      // adding interpolated value to data of cell
      daGausses[id]->SetTuple(j + offset, interps.data());
    }
  }
  return numGaussPoints;
}

operator=()

GaussCubature& GaussCubature::operator= ( const GaussCubature &  that )

delete

setArrayIDs()

void GaussCubature::setArrayIDs ( const std::vector< int > &  _arrayIDs )

inline

Definition at line 105 of file Cubature.H.

References cellType.

{ arrayIDs = _arrayIDs; }

writeGaussMesh()

void GaussCubature::writeGaussMesh

(

const char *

name

)

const

Definition at line 832 of file Cubature.C.

References gaussMesh.

{
  if (gaussMesh)
    writeVTFile<vtkXMLPolyDataWriter>(name, gaussMesh);
  else
  {
    std::cerr << "Gauss point mesh has not been constructed" << std::endl;
    exit(1);
  }
}

Member Data Documentation

arrayIDs

std::vector<int> GaussCubature::arrayIDs

private

Definition at line 145 of file Cubature.H.

Referenced by constructGaussMesh(), getGaussPointsAndDataAtCell(), integrateOverAllCells(), integrateOverCell(), and interpolateToGaussPoints().

dataSet

vtkSmartPointer<vtkDataSet> GaussCubature::dataSet

private

Definition at line 137 of file Cubature.H.

Referenced by constructGaussMesh(), GaussCubature(), getGaussPointsAndDataAtCell(), getOffset(), integrateOverAllCells(), integrateOverCell(), interpolateToGaussPoints(), and interpolateToGaussPointsAtCell().

dict

vtkQuadratureSchemeDefinition** GaussCubature::dict

private

Definition at line 143 of file Cubature.H.

Referenced by constructGaussMesh(), getGaussPointsAndDataAtCell(), integrateOverCell(), and interpolateToGaussPointsAtCell().

gaussMesh

vtkSmartPointer<vtkPolyData> GaussCubature::gaussMesh

private

Definition at line 141 of file Cubature.H.

Referenced by constructGaussMesh(), getGaussPointsAndDataAtCell(), integrateOverAllCells(), interpolateToGaussPoints(), and writeGaussMesh().

numComponents

std::vector<int> GaussCubature::numComponents

private

Definition at line 146 of file Cubature.H.

Referenced by getGaussPointsAndDataAtCell(), integrateOverAllCells(), and interpolateToGaussPoints().

numVolCells

int GaussCubature::numVolCells

private

Definition at line 139 of file Cubature.H.

Referenced by constructGaussMesh().

totalComponents

int GaussCubature::totalComponents

private

Definition at line 147 of file Cubature.H.

Referenced by getGaussPointsAndDataAtCell(), and interpolateToGaussPoints().

---

The documentation for this class was generated from the following files:

• Cubature.H
• Cubature.C