vtkHyperTreeGrid Paraview example in 2D

ParaView Support
1

Hi,

I was looking for some help to convert this Paraview Example to run in 2D. I’ve tried the obvious things and none of them seem to work. I’d like to have a 2D Quad that’s been refined; I’d like to leave the rest of the example as is, just make it work in 2D. Can anyone help? I’m using vtk (8.2.0) and paraview 5.8. I would also like to output the results, I’ve tried:

vtkSmartPointer<vtkXMLHyperTreeGridWriter>
    writer = vtkSmartPointer<vtkXMLHyperTreeGridWriter>::New();
writer->SetFileName("out.htg");
writer->SetInputData(hyperTreeGrid);
writer->Write();

But it’s always empty (cells and points) according to Paraview. To make this work in 2D, my question are:

vtkNew<vtkHyperTreeGrid> hyperTreeGrid;
hyperTreeGrid->Initialize();
int extent[6] = { 0, 1, 0, 1, 0, 1 }; // <===== Change something here, what?
hyperTreeGrid->SetGridExtent(extent);
hyperTreeGrid->SetDimension(3);  // <===== set to 2?
// hyperTreeGrid->SetOrientation(0);
hyperTreeGrid->SetBranchFactor(2);
// hyperTreeGrid->SetMaterialMaskIndex(nullptr);

vtkNew<vtkDoubleArray> xCoords;
xCoords->SetNumberOfValues(2);
xCoords->SetValue(0, 0);
xCoords->SetValue(1, 1.);
vtkNew<vtkDoubleArray> yCoords;
yCoords->SetNumberOfValues(2);
yCoords->SetValue(0, 0);
yCoords->SetValue(1, 1);

vtkNew<vtkDoubleArray> zCoords; 
zCoords->SetNumberOfValues(2);  // <===== Change to 1?
zCoords->SetValue(0, 0);
zCoords->SetValue(1, 1);  // <===== Remove?

hyperTreeGrid->SetXCoordinates(xCoords);
hyperTreeGrid->SetYCoordinates(yCoords);
hyperTreeGrid->SetZCoordinates(zCoords);

hyperTreeGrid->GenerateTrees();

Any thoughts? In particular how to write the mesh out once it runs in 2D.
Andy

2

Can anyone help? I tried the nightly release of paraview to see if it could see the content of the vtkXMLHyperTreeGridWriter but there is no difference. I tried rc3 of 9.0 and the changes to vtkHyperTreeGrid were so substantial that I could actually make my code run. Is there any chance somebody who uses vtkHyperTreeGrid could help; particularly for using it in 2D?
Thanks,
Andy

3

vtkHyperTreeGrid code changed a lot in the past 2 years. For example, vtkHyperTreeGrid used to inherit from vtkDataSet, it now inherits vtkDataObject, because its API and behavior are too different from vtkDataSet. And many of its API was also changed to defog a few ambiguities.

With the example you shared, you are generating a vtkHyperTreeGrid of one vtkHyperTree in 3D. I am not extremely familiar with this old API, but here are my advices to make it work in 2D: apply everything from your comments, and change extent to equal {0,1,0,1,0,0}.

Now, let’s go through a few features that might confuse you before you go in the wild.

If you want a grid of i x j x k vtkHyperTree, you need to set the coordinates to (i+1) x (j+1) x (k+1). This is because the coordinates that you give are the coordinates of the corners of the hypertrees. This is equivalent with considering the corners of the pixel of an image instead of the pixel itself, which could be mapped to the dual of a grid defined like that.

Now, there are no “points” or “cells” in a vtkHyperTreeGrid. There are trees, vertices, and leaves. Trees are to the vtkHyperTreeGrid what pixels are to an image. Vertices are the nodes of each tree, and leaves are the leaves of the trees. If you want to stick data to your vtkHyperTreeGrid, you can store it in a vtkDataArray which you feed to the vtkPointData obtained by calling hyperTreeGrid->GetPointData(). You can put multiple vtkDataArray in the same vtkPointData.

Adding data in your vtkHyperTreeGrid is not trivial. You have to run through it using what are called cursors, It is usually done by iterating over every trees, and by recursively refine the trees. You can find an example in ParaView master, in a plugin, at file Plugins/HyperTreeGridADR/HyperTreeGridFilters/vtkResampleToHyperTreeGrid.cxx. Look at the method vtkResampleToHyperTreeGrid::GenerateTrees and the recursive method vtkResampleToHyperTreeGrid::SubdivideLeaves. You should use the same kind of pattern.

FYI, you will see in current state of VTK or ParaView that data is not held anymore by a vtkPointData but by a vtkCellData, so don’t get confused by that. This change will not be in VTK 9.0 but in the next release.

4

If you want to learn more about the inner structure of vtkHyperTreeGrid and what families of cursors do what, you can check out this paper

5

Thank you very much for your reply and help. I had tried the extents as you suggest with two zeros in z but saw no difference. I wonder whether my 2D attempt been had correct all along and my problem is IO. I say this because if I take the paraview 3D FE example I referenced in my original post with no changes and add the following at the bottom (in 3D):

vtkNew<vtkXMLHyperTreeGridWriter> writer;
writer->SetFileName("new.htg");
writer->SetInputData(hyperTreeGrid);
writer->SetDataModeToAscii();
writer->Write();

The file is empty as far as paraview is concerned. Since this is the same effect I saw in trying to convert the example to 2D, I had just assumed the 2D was wrong. So to proceed I need to know how to write a hypertreegrid and read it in paraview. In 2/3D. Can you help?

I’ve reverted back to vtk 8.2 and paraview 5.8 for reference.

Thanks,
Andy

6

Can you share the file that you produced so I can take a peak at it?

7

Hi, Yep no problem. Code to produce file looks like:

... code from FE example

this->FillHTG(hyperTreeGrid);

vtkNew<vtkXMLHyperTreeGridWriter> writer;
writer->SetFileName("new.htg");
writer->SetInputData(hyperTreeGrid);
writer->SetDataModeToAscii();
writer->Write();

The file is attached. The above writer is just placed in the example after the FillHTG call. Looking at it on OS X on a stock download of paraview 5.8. Code ran with debug build of vtk 8.2 tag.

Thanks,
Andy

new.htg (207.9 KB)

8

Hi,
I wondered if the example was enough to figure out what was going on? Did you have a chance to look at this?
Thanks,
Andy

9

I quickly looked at it. I’ll have more insight on what is going on next week hopefully. But what I can say today is that vtkHyperTreeGrid was not a “stable” class in VTK 8.2. It was pretty much experimental. The core of the code changed a lot. It is getting more stable and consistent. There won’t be big API changes besides maybe method renaming. If you want to use this class in the long term, I would recommend using VTK 9.0 when it is released, or a release candidate in the meanwhile. The code is getting improved a lot (there were many design flaws in the former version), and it would be nice to not get stuck in the past on this data structure. You should take inspiration from vtkResampleToHyperTreeGrid plugin in ParaView 5.8 source code. This filter literally generates a hypertree grid from any input dataset.

10

@ajp I investigated, and it appears that the hypertree grid writer from VTK 8.2 is broken. It only worked on very specific situations. So you should update to VTK 9.0 and go from there. Things should work better.

11

Hi, I’m only just getting back to this now but thank you for your comments. I’ve recompiled against a clean build of vtk9. The example I posted in my first post does not work, as in really fails, headers have been removed etc. If you follow the post you can see all I really wanted to do was make the example I linked to work in 2D and write it out to file. The reason I’ve moved to vtk9 is based on your suggestion.

Can I confirm, none of those paraview examples are expected to work now and are deprecated?

The code posted in that example is exactly what I was looking to do, is there any way to reconstitute a similar example that works in vtk9 or is that type of thing not supported anymore?
Thanks,
Andy

12

The paraview example is indeed out of date, that will be updated.

There are 2 main ways you could create your htg.

  • You can do it “by hand”. Let me show an example in c++ (don’t forget to include relevant headers)
vtkNew<vtkHyperTreeGrid> htg;
/*
 * Initialize it like you did with x, y, and z coordinates
 */
vtkIdType treeOffset = 0;
for (vtkIdType i = 0; i < htg->GetCellDims()[0]; ++i)
{
  for (vtkIdType j = 0; j < htg->GetCellDims()[1]; ++j)
  {
     // We go through every trees in the grid. It is a 2D grid so there are 2 loop.
     // Let's initialize k = 0 for completeness.
    vtkIdType k = 0, treeId;
    // Get the tree index.
    htg->GetIndexFromLevelZeroCoordinates(treeId, i, j, k);
    // Create a cursor (there are many kinds, we can go through them later if you have questions)
    auto cursor = vtkSmartPointer<vtkHyperTreeGridNonOrientedCursor>::Take(
        htg->NewNonOrientedCursor(treeId, true /* create */);
    // Set tree indexing offset
    cursor->GetTree()->SetGlobalIndexStart(treeOffset);
    // You have to write your own procedure here, I assume you wrote the following
    // recursive function
    RecursivelySubdivideLeaves(cursor, paramsDrivingSubdivision);
    treeOffset += cursor->GetTree()->GetNumberOfVertices(); // we update the offset.
  }
}

Then, to subdivide a leave, you need to call cursor->SubdivideLeaf(). You have access of the global htg level index in the scalar fields with cursor->GetGlobalIndex(), and index relative to one tree with cursor->GetVertexId(). To navigate, you have cursor->GetNumberOfChildren() telling you how many children you have, cursor->ToChild(iChild) and cursor->ToParent() to walk inside the tree, and cursor->IsLeaf() to know when to stop.

  • You can also use vtkHyperTreeGridSource, which will create the structure of a htg given a descriptor. You can use it to create a htg pretty fast. The descriptor works as follow: . means “leaf”, R means subdivide, | means next level. An htg of 4 quadtrees like you showed in your first comment would be ..... If you want to subdivide the 4th tree, you would write ...R|..... If you want to subdivie the 2nd and the 4th tree, you would write .R.R|......... Finally, if you want to additionally subdivide the first child of the 4th tree, you would write .R.R|....R...|.....

Now, concerning writers, there might be a bug that went through the release. To make sure this works, you should call writer->SetDataSetMajorVersion(0) before your write. You don’t need to specify this parameter for the reader.

Did I answer your question or are there still some dark corners? Don’t hesitate to ask for clarification if something isn’t clear.

Thanks,
Yohann

13

@Yohann_Bearzi thank you very much for this. Your help has really allowed me to progress here! I’ve made a few advances including using vtkHyperTreeGridNonOrientedGeometryCursor as I wanted to obtain the bounds of the cell the cursor was pointing to. I wanted this so I could test whether an arbitrary point was contained within the cell the cursor was pointing to. It works, I assume it’s ok to do this from a performance perspective?

One thing I really do have a requirement to be able to do, and can’t find a cursor to allow me to do it, is depicted in the attached.

neighbours
neighbours1322×1306 8.32 KB
.

I would like to be able, for when the cursor points to the red cell (so may or may not be a leaf cell) to obtain/iterate over the red cell’s face neighbours, so those that have the outline in purple round them. 4 cell neighbours. Can you tell me how to do this please and also indicate if any specific method is more performant than another? I’d want to use the same approach in 3D also, just so you know.

If you are able to help me, I’d also like to:

  • Obtain the level for the cell the cursor points at. I think I know how to do this GetLevel, but wanted you to confirm, I presume this is zero numbered and is based on the number of times the root has been refined?
  • Obtain a unique, global id for each cell, that doesn’t change if any of the mesh is modified, only if that cell is removed.
  • How do I de-refine a cell, and all children below, can this be done? I’ve made some mistakes in refining the mesh where I didn’t need to, and currently I’m deleting the mesh and starting the whole process from scratch, it would be better to de-refine if possible.
  • How do I obtain the point locations for the 4/8 vertices that define the cell in 2D/3D respectively?

Many thanks again,
Andy

14

I’m happy you progressed! I’ll answer the questions point by point:

  • Yes
  • GetGlobalIndex() in the cursor is unique.
  • You cannot de-refine a cell in the current implementation. However you can mask them. The vtkBitArray used for the mask is indexed by GetGlobalIndex in the cursor. When you mask a node, it is discarded from display, as well as all of its subtree. We have this mechanism because since htg have implicit indexing that is constructed as you refine, you would break everything is you could remove a cell. With masking you can do the same effect at the cost of one bit.
  • You should be able to catch the points coordinates using the bounds available in the geometry cursor. I’d say that the half-sum in each dimension gives you the center :slight_smile:
15

Ok thanks for the reply. Could I ask for you thoughts/advice regarding the question about face neighbours?

I would like to be able, for when the cursor points to the red cell (so may or may not be a leaf cell) to obtain/iterate over the red cell’s face neighbours, so those that have the outline in purple round them. 4 cell neighbours. Can you tell me how to do this please and also indicate if any specific method is more performant than another? I’d want to use the same approach in 3D also, just so you know

Thanks,
Andy

16

Didn’t catch this question. You need to use a super cursor for this kind of work. vtkHyperTreeGridNonOrientedVonNeumannSuperCursor gives you a 4-neighborhood in 2D. Replace VonNeumann by Moore and you’ll have a 8-neighborhood. You can get a neighborhood cursor by calling GetOrientedGeometryCursor(iCursor). Be careful however, in 2D with a VonNeumann cursor, index 2 returns yourself, the center cursor. Don’t try to get non oriented geometry cursors using this method as they do not work (we might deprecate these methods soon). An oriented cursor does not have access to ToParent()

17

Last thing: look at the available methods in supercursors, in particular in vtkHyperTreeGridNonOrientedSuperCursor, as you might not need to create a cursor at your neighbor’s to access to most information you might need.

18

Ok, I’ll take a look at all of this.
Thanks
Andy

19

@Yohann_Bearzi I was wondering if you could help again, I’m pretty stuck. I’m trying to find all of the neighbour cells (in green) that border the red cell. Please see below:

neigh
neigh2634×1370 21.4 KB

My plan was for example to look to the west neighbours of the red cell, traverse the tree until I get to the leafs, but use the fact that the east neighbour of those leaf cells must match my red cell for them to be touching. Sort of a bit lost. This is what I’ve got so far, but I can’t finish the traversal bit:

constexpr int SOUTH = 1;
constexpr int WEST = 3;
constexpr int EAST = 5;
constexpr int NORTH = 7;

void NeighbourProcess(vtkHyperTreeGridOrientedGeometryCursor *cursor, const int & /*direction*/, std::vector<std::array<double, 3>> &neighbourCenters)
{
    if (cursor->HasTree())
    {
        if (cursor->IsLeaf())
        {
            std::array<double, 3> pt;
            cursor->GetPoint(pt.data());
            neighbourCenters.push_back(pt);
        }
        else
        {
            //????
        }
    }
}

void Compute(vtkHyperTreeGridNonOrientedMooreSuperCursor *superCursor)
{
    if (superCursor->IsLeaf())
    {
        // https://gitlab.kitware.com/todoooo/vtk/blob/cf111cf65e200acabfb87188867a078ec0da6035/Filters/HyperTree/vtkHyperTreeGridToDualGrid.cxx
        std::vector<std::array<double, 3>> neighbourCenters;
        vtkSmartPointer<vtkHyperTreeGridOrientedGeometryCursor> cursorS = superCursor->GetOrientedGeometryCursor(SOUTH);
        NeighbourProcess(cursorS, NORTH, neighbourCenters);
        vtkSmartPointer<vtkHyperTreeGridOrientedGeometryCursor> cursorW = superCursor->GetOrientedGeometryCursor(WEST);
        NeighbourProcess(cursorW, EAST, neighbourCenters);
        vtkSmartPointer<vtkHyperTreeGridOrientedGeometryCursor> cursorE = superCursor->GetOrientedGeometryCursor(EAST);
        NeighbourProcess(cursorE, WEST, neighbourCenters);
        vtkSmartPointer<vtkHyperTreeGridOrientedGeometryCursor> cursorN = superCursor->GetOrientedGeometryCursor(NORTH);
        NeighbourProcess(cursorN, SOUTH, neighbourCenters);
    }
    else
    {
        const int numberOfChildren = static_cast<int>(superCursor->GetNumberOfChildren());
        for (int ichild = 0; ichild < numberOfChildren; ++ichild)
        {
            superCursor->ToChild(static_cast<unsigned char>(ichild));
            Compute(superCursor);
            superCursor->ToParent();
        }
    }
}

void Process(vtkHyperTreeGrid *htg)
{
    vtkHyperTreeGrid::vtkHyperTreeGridIterator it;
    htg->InitializeTreeIterator(it);
    vtkIdType treeId;
    while (it.GetNextTree(treeId))
    {
        vtkNew<vtkHyperTreeGridNonOrientedMooreSuperCursor> superCursor;
        superCursor->Initialize(htg, treeId);
        Compute(superCursor);
    }
}

For this bit of test code I was initially planning to store the centre of the neighbours I found. I can then output that for a target cell and look at the points in paraview to confirm it’s correct.

Eventually, what I actually want to find are the nodes that “define” the red cell, that is not just the base 4 corner points, but all the nodes along the refined edges. To be able to do this the way I planned I need all of the green neighbours of the red cell. I could then put the 4 corner points of each green neighbour into a pointlocator, and store the index of the point in a multiset. Those points with a count == 2 in addition to the 4 base corner points of the red cell would result in a list of points defining the red cell. I can’t really think of another way to do it. I want to output that cell as a polygon.

I don’t suppose there is any easy function to give me the list of refined nodes that make up the refined edges of the original 4 base edges of the red cell?

Many thanks,
Andy

20

To do that, you’ll want to use supercursors (the Von Neumann one for your particular case) (see this paper for reference). You can call GetOrientedGeometryCursor(iChild) to create a vtkOrientedGeometryCursor (don’t use the non oriented version). From this cursor, if you were at the north cursor, you go south and extract the leafs interfacing the south. and the same thing goes for every directions.