1
2
3
4
5
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
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
|
// SPDX-FileCopyrightText: Copyright (c) Kitware Inc.
// SPDX-License-Identifier: BSD-3-Clause
#include "FEDataStructures.h"
#include <cstdlib>
#include <iostream>
#include <mpi.h>
#include <vtkCPInputDataDescription.h>
#include <vtkCellData.h>
#include <vtkCellType.h>
#include <vtkDoubleArray.h>
#include <vtkExtentTranslator.h>
#include <vtkImageData.h>
#include <vtkNew.h>
#include <vtkPointData.h>
#include <vtkUnsignedCharArray.h>
#include <vtkUnstructuredGrid.h>
Grid::Grid(const unsigned int numberOfPoints[3], bool outputImageData, int numberOfGhostLevels)
{
if (numberOfPoints[0] == 0 || numberOfPoints[1] == 0 || numberOfPoints[2] == 0)
{
std::cerr << "Must have a non-zero amount of points in each direction.\n";
}
int wholeExtent[6];
for (int i = 0; i < 3; i++)
{
wholeExtent[2 * i] = 0;
wholeExtent[2 * i + 1] = numberOfPoints[i] - 1;
}
// we use the ExtentTranslator to compute the partitioning
int mpiSize = 1;
int mpiRank = 0;
MPI_Comm_rank(MPI_COMM_WORLD, &mpiRank);
MPI_Comm_size(MPI_COMM_WORLD, &mpiSize);
vtkNew<vtkExtentTranslator> extentTranslator;
int extent[6];
extentTranslator->PieceToExtentThreadSafe(
mpiRank, mpiSize, numberOfGhostLevels, wholeExtent, extent, 3, 0);
// create the grid
if (outputImageData)
{
this->CreateImageData(extent);
}
else
{
this->CreateUnstructuredGrid(extent);
}
// set ghost level information...
vtkNew<vtkUnsignedCharArray> ghostCells;
ghostCells->SetNumberOfTuples(this->VTKGrid->GetNumberOfCells());
ghostCells->SetName(vtkDataSetAttributes::GhostArrayName());
ghostCells->Fill(0);
this->VTKGrid->GetCellData()->AddArray(ghostCells);
if (numberOfGhostLevels || mpiSize > 1)
{
// get the extent if we didn't have ghost levels so that we can compare
// which cells are ghost cells
int nonGhostLevelExtent[6];
extentTranslator->PieceToExtentThreadSafe(
mpiRank, mpiSize, 0, wholeExtent, nonGhostLevelExtent, 3, 0);
vtkIdType counter = 0;
for (int k = extent[4]; k < extent[5]; k++)
{
bool zGhost = k < nonGhostLevelExtent[4] || k >= nonGhostLevelExtent[5];
for (int j = extent[2]; j < extent[3]; j++)
{
bool yGhost = j < nonGhostLevelExtent[2] || j >= nonGhostLevelExtent[3];
for (int i = extent[0]; i < extent[1]; i++)
{
bool xGhost = i < nonGhostLevelExtent[0] || i >= nonGhostLevelExtent[1];
if (xGhost || yGhost || zGhost)
{
unsigned char value = vtkDataSetAttributes::DUPLICATECELL;
ghostCells->SetTypedTuple(counter, &value);
}
counter++;
}
}
}
}
}
vtkDataSet* Grid::GetVTKGrid()
{
return this->VTKGrid;
}
void Grid::CreateImageData(int extent[6])
{
vtkNew<vtkImageData> imageData;
imageData->SetExtent(extent);
this->VTKGrid = imageData;
}
void Grid::CreateUnstructuredGrid(int extent[6])
{
// create the points -- slowest in the x and fastest in the z directions
double coord[3] = { 0, 0, 0 };
vtkNew<vtkPoints> points;
for (int i = extent[0]; i <= extent[1]; i++)
{
coord[0] = static_cast<double>(i);
for (int j = extent[2]; j <= extent[3]; j++)
{
coord[1] = static_cast<double>(j);
for (int k = extent[4]; k <= extent[5]; k++)
{
coord[2] = static_cast<double>(k);
points->InsertNextPoint(coord);
}
}
}
vtkNew<vtkUnstructuredGrid> unstructuredGrid;
unstructuredGrid->Initialize();
unstructuredGrid->SetPoints(points);
// create the hex cells
vtkIdType cellPoints[8];
int numPoints[3] = { extent[1] - extent[0] + 1, extent[3] - extent[2] + 1,
extent[5] - extent[4] + 1 };
unstructuredGrid->Allocate((numPoints[0] - 1) * (numPoints[1] - 1) * (numPoints[2] - 1));
for (int k = 0; k < numPoints[2] - 1; k++)
{
for (int j = 0; j < numPoints[1] - 1; j++)
{
for (int i = 0; i < numPoints[0] - 1; i++)
{
cellPoints[0] = i * numPoints[1] * numPoints[2] + j * numPoints[2] + k;
cellPoints[1] = (i + 1) * numPoints[1] * numPoints[2] + j * numPoints[2] + k;
cellPoints[2] = (i + 1) * numPoints[1] * numPoints[2] + (j + 1) * numPoints[2] + k;
cellPoints[3] = i * numPoints[1] * numPoints[2] + (j + 1) * numPoints[2] + k;
cellPoints[4] = i * numPoints[1] * numPoints[2] + j * numPoints[2] + k + 1;
cellPoints[5] = (i + 1) * numPoints[1] * numPoints[2] + j * numPoints[2] + k + 1;
cellPoints[6] = (i + 1) * numPoints[1] * numPoints[2] + (j + 1) * numPoints[2] + k + 1;
cellPoints[7] = i * numPoints[1] * numPoints[2] + (j + 1) * numPoints[2] + k + 1;
unstructuredGrid->InsertNextCell(VTK_HEXAHEDRON, 8, cellPoints);
}
}
}
this->VTKGrid = unstructuredGrid;
}
void Grid::UpdateField(double time, vtkCPInputDataDescription* inputDataDescription)
{
if (inputDataDescription->IsFieldNeeded("Scalar", vtkDataObject::POINT))
{
vtkDoubleArray* scalar =
vtkDoubleArray::FastDownCast(this->VTKGrid->GetPointData()->GetArray("Scalar"));
if (scalar == nullptr)
{
scalar = vtkDoubleArray::New();
scalar->SetNumberOfTuples(this->VTKGrid->GetNumberOfPoints());
scalar->SetName("Scalar");
this->VTKGrid->GetPointData()->AddArray(scalar);
scalar->Delete(); // ok since VTKGrid is keeping a reference to this array
}
vtkIdType numPoints = this->VTKGrid->GetNumberOfPoints();
for (vtkIdType pt = 0; pt < numPoints; pt++)
{
double coord[3];
this->VTKGrid->GetPoint(pt, coord);
double value = coord[1] * time;
scalar->SetTypedTuple(pt, &value);
}
}
}
|
