Memory layout of blocks

[CsatiZoltan]

Hi,

My goal is to construct a vtkUnstructuredGrid object from a libMeshb mesh without making copy. For that, I need to understand the memory layout of libMeshb.

I use GmfGetBlock to get the data in one chunk. For now, I consider only the vertices, but I will also have to fetch the elements (currently, we are only interested in low order Lagrange elements). Based on Section 3.8 of the documentation, we have two options:

1. Get the vertices by coordinates, in separate arrays. This is the scalar approach, resulting in contiguous arrays x, y and z (in 3D). I am not sure that I can pass them separately to vtkPoints.
2. Get the vertices in a single variable, called vector addresses. Looking at the documentation and the example code in
   test_libmeshb_block.c, the vertices are stored in an array of pointers. My issue with this approach is that the whole table is not contiguous in memory, so I cannot expose it as a buffer to VTK. Is there a reason you do not store them in a single array? What about storing them as [x1, y1, z1, x2, y2, z2, ..., xN, yN, zN] or [x1, x2, ..., xN, y1, y2, ..., yN, z1, z2, ..., zN] ? The element-vertex connectivities could also be stored in a single array, given that the mesh format stores each element type separately.

By the way, what is the role of the reference as the last entry a vertex or an element @LoicMarechal ?

[LoicMarechal]

My goal is to construct a vtkUnstructuredGrid object from a libMeshb mesh without making copy. For that, I need to understand the memory layout of libMeshb. I use GmfGetBlock to get the data in one chunk. For now, I consider only the vertices, but I will also have to fetch the elements (currently, we are only interested in low order Lagrange elements). Based on Section 3.8 of the documentation, we have two options: • Get the vertices by coordinates, in separate arrays. This is the scalar approach, resulting in contiguous arrays x, yand z (in 3D). I am not sure that I can pass them separately to vtkPoints. • Get the vertices in a single variable, called vector addresses. Looking at the documentation and the example code in test_libmeshb_block.c, the vertices are stored in an array of pointers. My issue with this approach is that the whole table is not contiguous in memory, so I cannot expose it as a buffer to VTK. Is there a reason you do not store them in a single array? What about storing them as [x1, y1, z1, x2, y2, z2, ..., xN, yN, zN] or [x1, x2, ..., xN, y1, y2, ..., yN, z1, z2, ..., zN] ? The element-vertex connectivities could also be stored in a single array, given that the mesh format stores each element type separately. By the way, what is the role of the reference as the last entry a vertex or an element?

Hi Zoltan GmfGetBlock() does have any internal layout but uses the memory layout you provide. If you have a single COORD table that stores x1, y1, z1, x2, y2, z2 , ... You can call: GmfGetBlock(InpMsh, GmfVertices, 1, NmbVer, 0, NULL, NULL, GmfFloatVec, 3, &COORD[ 1 * 3 ], &COORD[ NmbVer * 3 ], GmfInt, &RefTab[1], &RefTab[ NmbVer ] ); It also works with coordinates stored inside some kind of structures or objects : // Define a vertex type typedef struct { int MyData; float coords[3[; int ref; ... }MyVertex // Define a table of 1000 vertices MyVertex Vertices[1000]; // Read the vertices from file GmfGetBlock(InpMsh, GmfVertices, 1, NmbVer, 0, NULL, NULL, GmfFloatVec, 3, &Vertices[1].coords[0], &Vertices[ NmbVer ].coords[0], GmfInt, &Vertices[1].ref, &Vertices[ NmbVer ].ref); Even if the sets of three coordinates are not contiguous in memory, the libMeshb will compute the memory stride and jump through memory to the right address. The last integer is a vertex arbitrary attribute you can use freely to store some physical attribute or a patch reference, for example. Hope it helps, Loïc

[CsatiZoltan]

Thank you for the explanations.

Your first method (COORD array) is a proper choice for me: I store the x, y, (z) coordinates in an std::vector<double>, which works fine because the underlying buffer is guaranteed by the C++ standard to be contiguous. My problem was that I forgot changing GmfFloatVec to GmfDoubleVec when pasting the example code.

Even if the sets of three coordinates are not contiguous in memory, the libMeshb will compute the memory stride and jump through memory to the right address.

libMeshb may compute the strides, but external libraries are often not flexible enough to allow the user to provide the stride, so they usually require a contiguous chunk of memory. This is why I insisted on contiguity.

Since then, I have another question. Why do you start the reading from index 1? In my experience, it reads a row full of zeros. So instead of allocating storage for (NmbVer + 1)*3 elements and calling

&COORD[ 1 * 3 ], &COORD[ NmbVer * 3 ]

I allocate storage for NmbVer*3 elements and call

&COORD[0], &COORD[ NmbVer * 3 - 3 ]

The same reasoning goes for the scalar address, e.g. for the references. I.e.

GmfGetBlock(idx, GmfVertices, 1, nVertex, 0, nullptr, nullptr,
            GmfDoubleVec, 3, &COORD[0],  &COORD[3*NmbVer-3],
            GmfInt,          &RefTab[0], &RefTab[NmbVer-1]);

It works fine for me.
For a simple example, consider the following mesh file:

MeshVersionFormatted 2

Dimension 3

Vertices
3
0.0500000007451 0.0046905297786 0.00167893001344 0
0.0500000007451 0.00427286000922 0.00256174011156 0
0.0500000007451 0.00341648003086 0.00201285001822 1

End

With your method, I get

COORDS=[
0
0
0
0.050000000745099998
0.0046905297786000002
0.0016789300134399999
0.050000000745099998
0.0042728600092199996
0.0025617401115599998
0.050000000745099998
0.0034164800308599998
0.0020128500182199998]

With my way, I obtain

COORDS=[
0.050000000745099998
0.0046905297786000002
0.0016789300134399999
0.050000000745099998
0.0042728600092199996
0.0025617401115599998
0.050000000745099998
0.0034164800308599998
0.0020128500182199998]

[LoicMarechal]

Since then, I have another question. Why do you start the reading from index 1? In my experience, it reads a row full of zeros. So instead of allocating storage for (NmbVer + 1)*3 elements and calling

It was just for the sake of example. Some people like to start from 1 to N, some from 0 to N-1. You can send anything you like to the libMeshb as long as there is the right number of vertices. You can even use the libMeshb in a parallel code, with each process reading the same file at the same time but accessing different blocks of data like that: GmfGetBlock(idx, GmfVertices, BeginIndex, EndIndex, 0, nullptr, nullptr, GmfDoubleVec, 3, &COORD[ BeginIndex * 3 ], &COORD[ EndIndex * 3 ], GmfInt, &RefTab[ BeginIndex ], &RefTab[ EndIndex ]); Loïc

[CsatiZoltan]

You can even use the libMeshb in a parallel code, with each process reading the same file at the same time but accessing different blocks of data like that

Would it work with distributed parallelism? It would be very interesting for large meshes: each MPI process would read a chunk and pass that part of the mesh to the ParaView server (which also supports distributed parallelism). The elements could be partitioned, but the vertices should be present on every process, I guess.