VTK和ParaView中引入了显式结构化网格表达地质网格
Introducing Explicit Structured Grids in VTK and ParaView - Kitware Blog
1、简介
新版本的vtk引入了适用于油藏角点网格模型的显式结构化网格，即在原来非结构化网格vtkUnstructuredGrid和结构化网格vtkStructuredGrid的基础上添加了处于两者之间的显式结构化网格vtkExplicitStructuredGrid。它表示有拓扑结构的数据集（所有单元格均为六面体，并沿i，j和k方向进行结构化），但是可以显式地定义几何特点。该数据集可以表示地质或储层网格，例如在任何方向都可能存在断层的地方。
图1Johansen formation 地质模型公开数据用这个方式显示的效果https://www.sintef.no/projectweb/matmora/downloads/johansen/
2、数据结构
显式结构化网格是VTK中的一种新数据结构，可共享结构化和非结构化网格的属性。 这个新的数据集非常适合于表示在油气领域众所周知的储层网格。
显式结构化网格的数据结构使用结构化坐标系（i，j，k）确定单元邻域，这也支持单元消隐。 同时，每个单元的几何形状由其顶点明确定义（请参见图2）。 这种定义单元的方式不仅可以定义曲线网格，而且可以定义具有更复杂拓扑的网格。 在明确的结构化网格中，单元可以相邻，但不需要连接。 这些几何形状上的不连续性可以以与非结构化网格相同的效率对单元之间的大小裂缝进行建模，而无需增加单元分辨率。 但是，与非结构化网格相比，显式结构化网格具有其底层规则结构的优势，可以在数据并行上下文中获得更好的性能（请参见图3）。
图2 三种网格的比较
图3 三种网格优点缺点对比表
3、过滤器
目前，这种新的数据结构有优化过的截取器和集中转化器，以及用于显式网格的算法的类，转换器使用户能够轻松地将其现有数据转换为显式的结构化网格，从而使用户可以更轻松地访问新数据结构。
这些转换器包括：
To / From UnstructuredGrid；
From ImageData;
To PolyData（即提取曲面网格）;
显式结构化网格生成器源。
最后，一个新的数据源也被添加到Paraview中，以快速生成不同的显式结构化网格储层数据（参见图4）。
图4 显示结构网格生成源的5种油藏模型生成模式
4、数据格式
由于包含了转换过滤器以及与非结构化网格的数据结构相似性，因此不需要特定于显式结构化网格的数据格式。 相反，显式结构化网格将其数据保存为非结构化网格，而不会丢失任何信息。
5、实例测试
利用vtkExplicitStructutredGrid显示模型需要自定义每个网格八个点的坐标，每个点的属性值，每个六面体网格的定义，需要用到vtkPoints, vtkFloatArray, vtkHexahedron, vtkCellArray,具体代码及显示效果如下。
#include <iostream>
#include <random>
#include <armadillo>
#include <vtk-9.0/vtkType.h>
#include <vtk-9.0/vtkPoints.h>
#include <vtk-9.0/vtkPointData.h>
#include <vtk-9.0/vtkFloatArray.h>
#include <vtk-9.0/vtkHexahedron.h>
#include <vtk-9.0/vtkCellArray.h>
#include <vtk-9.0/vtkExplicitStructuredGrid.h>
#include <vtk-9.0/vtkDataSetMapper.h>
#include <vtk-9.0/vtkActor.h>
#include <vtk-9.0/vtkProperty.h>
#include <vtk-9.0/vtkScalarBarWidget.h>
#include <vtk-9.0/vtkScalarBarActor.h>
#include <vtk-9.0/vtkLookupTable.h>
#include <vtk-9.0/vtkRenderer.h>
#include <vtk-9.0/vtkRenderWindow.h>
#include <vtk-9.0/vtkRenderWindowInteractor.h>
#include <vtk-9.0/vtkAutoInit.h>
using namespace arma;
using namespace std;
struct Point3D
{
float xyz[3];
Point3D() {};
Point3D(float x, float y, float z) {
this->xyz[0] = x;
this->xyz[1] = y;
this->xyz[2] = z;
}
void setValue(Point3D p)
{
xyz[0] = p.xyz[0];
xyz[1] = p.xyz[1];
xyz[2] = p.xyz[2];
}
};
struct Hexahedron {
Point3D _p1, _p2, _p3, _p4, _p5, _p6, _p7, _p8;
Hexahedron() {};
Hexahedron(Point3D p1, Point3D p2, Point3D p3, Point3D p4,
Point3D p5, Point3D p6, Point3D p7, Point3D p8)
{
this->_p1.setValue(p1);
this->_p2.setValue(p2);
this->_p3.setValue(p3);
this->_p4.setValue(p4);
this->_p5.setValue(p5);
this->_p6.setValue(p6);
this->_p7.setValue(p7);
this->_p8.setValue(p8);
}
Hexahedron(Point3D p0, float dx, float dy, float dz) {
this->_p1.setValue(p0);
this->_p2.setValue(Point3D(p0.xyz[0] + dx, p0.xyz[1], p0.xyz[2]));
this->_p3.setValue(Point3D(p0.xyz[0] + dx, p0.xyz[1] + dy, p0.xyz[2]));
this->_p4.setValue(Point3D(p0.xyz[0], p0.xyz[1] + dy, p0.xyz[2]));
this->_p5.setValue(Point3D(p0.xyz[0], p0.xyz[1], p0.xyz[2] + dz));
this->_p6.setValue(Point3D(p0.xyz[0] + dx, p0.xyz[1], p0.xyz[2] + dz));
this->_p7.setValue(Point3D(p0.xyz[0] + dx, p0.xyz[1] + dy, p0.xyz[2] + dz));
this->_p8.setValue(Point3D(p0.xyz[0], p0.xyz[1] + dy, p0.xyz[2] + dz));
}
};
VTK_MODULE_INIT(vtkRenderingOpenGL2);
VTK_MODULE_INIT(vtkInteractionStyle);
VTK_MODULE_INIT(vtkRenderingFreeType);
int main()
{
std::cout << "Hello World!\n";
default_random_engine e;
uniform_real_distribution<float> rand(0, 1);
int ni = 100;
int nj = 100;
int nk = 50;
float lengthX = 10000;
float lengthY = 10000;
float lengthZ = 1000;
float dx = lengthX / ni;
float dy = lengthY / nj;
float dz = lengthZ / nk;
fmat bot = fmat(nj, ni, fill::zeros);
field<Hexahedron> gridPoints = field<Hexahedron>(nj, ni, nk);
vtkSmartPointer<vtkPoints> points = vtkPoints::New();
vtkSmartPointer<vtkFloatArray> scalars = vtkFloatArray::New();
//create the vtkExplicitStructuredGrid
vtkSmartPointer<vtkExplicitStructuredGrid> grid = vtkExplicitStructuredGrid::New();
grid->SetDimensions(ni, nj, nk);
vtkSmartPointer<vtkCellArray> cellArray = vtkCellArray::New();
vtkIdType seq = 0;
for (int k = 0; k < nk; k++) {
for (int j = 0; j < nj; j++) {
for (int i = 0; i < ni; i++) {
Point3D p0 = Point3D(i * dx, j * dy, bot(j, i) + dz * k);
gridPoints(j, i, k) = Hexahedron(p0, dx, dy, dz);
float value = rand(e);
points->InsertNextPoint(gridPoints(j, i, k)._p1.xyz);
scalars->InsertNextTuple1(value);
points->InsertNextPoint(gridPoints(j, i, k)._p2.xyz);
scalars->InsertNextTuple1(value);
points->InsertNextPoint(gridPoints(j, i, k)._p3.xyz);
scalars->InsertNextTuple1(value);
points->InsertNextPoint(gridPoints(j, i, k)._p4.xyz);
scalars->InsertNextTuple1(value);
points->InsertNextPoint(gridPoints(j, i, k)._p5.xyz);
scalars->InsertNextTuple1(value);
points->InsertNextPoint(gridPoints(j, i, k)._p6.xyz);
scalars->InsertNextTuple1(value);
points->InsertNextPoint(gridPoints(j, i, k)._p7.xyz);
scalars->InsertNextTuple1(value);
points->InsertNextPoint(gridPoints(j, i, k)._p8.xyz);
scalars->InsertNextTuple1(value);
seq += 8;
vtkSmartPointer < vtkHexahedron > cell = vtkSmartPointer<vtkHexahedron>::New();
cell->GetPointIds()->SetId(0, seq - 8);
cell->GetPointIds()->SetId(1, seq - 7);
cell->GetPointIds()->SetId(2, seq - 6);
cell->GetPointIds()->SetId(3, seq - 5);
cell->GetPointIds()->SetId(4, seq - 4);
cell->GetPointIds()->SetId(5, seq - 3);
cell->GetPointIds()->SetId(6, seq - 2);
cell->GetPointIds()->SetId(7, seq - 1);
cellArray->InsertNextCell(cell);
}
}
}
grid->SetPoints(points);
grid->GetPointData()->SetScalars(scalars);
grid->SetCells(cellArray);
//Create a mapper
vtkSmartPointer<vtkDataSetMapper> mapper = vtkDataSetMapper::New();
mapper->SetInputData(grid);
double* range;
range = scalars->GetRange();
cout << "cellColor:min=" << range[0] << " max=" << range[1] << endl;
mapper->SetScalarRange(range[0], range[1]);
vtkSmartPointer<vtkRenderWindowInteractor> renWinIn = vtkRenderWindowInteractor::New();
//create scalarbar
vtkSmartPointer<vtkLookupTable> lut = vtkSmartPointer<vtkLookupTable>::New();
lut->SetNumberOfColors(256);
lut->SetHueRange(0.667, 0);
//lut->SetRange(range[0], range[1]); //value below 0.01
lut->SetRange(0, 1); //value below 0.01
lut->Build();
mapper->SetLookupTable(lut);
vtkSmartPointer< vtkScalarBarActor > scalarBarActor = vtkSmartPointer< vtkScalarBarActor >::New();
scalarBarActor->SetOrientationToHorizontal();
scalarBarActor->SetLookupTable(lut);
vtkSmartPointer< vtkScalarBarWidget > scalarBarWidget = vtkSmartPointer< vtkScalarBarWidget >::New();
scalarBarWidget->SetInteractor(renWinIn);
scalarBarWidget->SetScalarBarActor(scalarBarActor);
scalarBarWidget->On();
//Create the actor
vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
actor->SetMapper(mapper);
actor->GetProperty()->SetDiffuseColor(1, 1, 0);
vtkSmartPointer<vtkRenderer> render = vtkRenderer::New();
vtkSmartPointer<vtkRenderWindow> renWin = vtkRenderWindow::New();
//Set the interactive style
renWinIn->SetRenderWindow(renWin);
renWin->AddRenderer(render);
render->SetBackground(0, 0, 0);
render->AddActor(actor);
render->AddActor(scalarBarActor);
render->ResetCamera();
renWin->Render();
renWinIn->Start();
}
