vtkHyperTreeGrid Paraview example in 2D

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)

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

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.

@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.

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

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

@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.

.

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

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:

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

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()

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.

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

@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:


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

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.

Thanks. That’s pretty consistent with the code I posted above. A question is:

From this cursor, if you were at the north cursor, you go south and extract the leafs interfacing the south

I posted code trying to do that. How do I know that a random child of the neighbour cell, borders the cell I’m interested in? So I know if I’m on the north face I want to compare to the south cell neighbours, how do I get the array of south cell neighbours? I can’t figure that out.

Is there an easier way just to extract all the points the define the perimeter of the red cell? That’s all I’m really trying to do

Right! Now I see your problem…

As an aside:

Why is Moore used Here for the 2D case and not Von Neumann? This is the reference I’ve been using for applicable examples. That’s why my code above uses Moore

You would want to use a Moore cursor if you wanted to turn the bottom left of the screen you shared green. Moore cursor is a 8-neighborhood in 2D (or the unit-disk using the infinity norm), Von Neumann is 4-neighborhood (or the unit-disk using the 0-norm).

Alright, you cannot do the perfect traversal in the current state of the htg library. What you will need to do as a fix in the meanwhile is clone your cursor, calling the Clone method before you go to the children. We’ll come with a better design in the future.

Ideally you would need to create a non oriented geometry cursor. It doesn’t work (if you look at the code, you’ll see an assert(false) which was put here to prevent people from using it) because it would be very expensive to create an actual non oriented cursor having the ability to go to the root of the tree. Ideally, I think we would need to have a cursor that cannot go to the true root, but traverse the subtree freely without the possibility to go to the parent of the neighbor.

We’ll address that, hopefully in the next release.

Ok thank again. I’m not sure I’m going to be able to get this working, I can’t really follow how to collect a list of neighbour cells that touch the perimeter of the cell that interests me; red cell in this case.

Can I just check: Is there an easier way just to extract all the points the define the perimeter of the red cell? That’s all I’m really trying to do, ignoring the code and approach I had posted, I’m open to any suggestions. Upon knowing the perimeter nodes I should be able to turn the red cell into a polygon.