Handling pre-built structured grids

[banesullivan]

the structured methods in this library (to my knowledge) only support 1-D arrays representing the axial coordinates of a rectilinear grid (the inputs to the numpy.meshgrid) for structured points. This is okay... but in cases of curvilinear structured grids where the user has their own tooling for building those grids, we may want to pass prebuilt coordinate arrays.

For example, take this curvilinear structured grid: curvi_mesh.vts.zip

import pyvista as pv

mesh = pv.read("curvi_mesh.vts")

mesh.plot(show_edges=True, show_grid=True)

This type of grid cannot be made with NumPy's meshgrid function as it cannot be represented as 3 1D coordinate arrays. This is a set of structured coordinates just like what numpy.meshgrid yields. These structured coordinates can be accessed by the x, y, and z attributes of a PyVista StructuredGrid object.

Note that this is something that is really common outside of PyVista too - take discretize's CurvilinearMesh class for instance.

>>> mesh.x.shape, mesh.y.shape, mesh.z.shape
((17, 17, 5), (17, 17, 5), (17, 17, 5))

So how can we pass these kinds of pre-built structured points to the field classes such that we can leverage any performance/optimization for structured grids and plotting routines?

Just to drive home that this is a pretty common thing, even matplotlib can handle these types of structured coordinates:

plt.pcolormesh(mesh.x[:,:,0], 
               mesh.y[:,:,0], 
               np.random.rand(*mesh.dimensions[0:2]),
               color="k")
plt.axis("image")
plt.show()

---

So basically, how might I use that structured gird with GSTools without just sending the points as an unstructured pointset. e.g.

from gstools import SRF, Gaussian

model = Gaussian(dim=2, var=0.2, len_scale=0.05)
srf = SRF(model, seed=20170519)

field = srf.structured([mesh.x[:,:,0], 
                        mesh.y[:,:,0]])

>>> field.shape
(289, 289)

This is the wrong shape as it should be (17, 17) but it instead flattened the x and y arrays that I passed and used all of those coordinates... (17*17=289)

srf.plot()

This is not the plot that I would expect.

[banesullivan]

Also, is there no advantage for the kriging when using structured pointsets? It looks like structured grids are cast to unstructured when solving the kriging problem...what exactly is happening here?

GSTools/gstools/krige/ordinary.py

Lines 88 to 94 in 5326163

	if mesh_type == "structured":
	mesh_type_changed = True
	mesh_type_old = mesh_type
	mesh_type = "unstructured"
	x, y, z, axis_lens = reshape_axis_from_struct_to_unstruct(
	self.model.dim, x, y, z
	)

[banesullivan]

Just FYI... casting to an unstructured point set works fine:

field = srf.unstructured([mesh.x[:,:,0].ravel(), 
                          mesh.y[:,:,0].ravel()])
srf.plot()
plt.axis("image")
plt.show()

So really, my question revolves around whether using "structured" points has any perfromance advantages and if so, how do we use custom structured pointsets? If there's no advantage, then all of this is irrelavant and we can just send these kinds of curvilinear grids as unstructured pointsets.

[MuellerSeb]

Structured fields give some advantages in the SRF class, since the cython generator has an optimized routine for that case.

I know, that coming from VTK, the term "structured" means a lot there. But we only support rectilinear grids for our structured meshes. In general it is best to use unstructured grids as you supposed, that is why this is the standard value for mesh_type in the __call__ routine.

I think it would be nice, if we could generate fields directly on pyvista meshes like you have proposed in #59 where these meshes are always recasted to unstructured meshes. Therefore it would also be nice to have a better handling of different generated fields in the class.

[banesullivan]

But we only support rectilinear grids for our structured meshes. In general it is best to use unstructured grids as you supposed, that is why this is the standard value for mesh_type in the call routine.

Can I recommend changing the term "structured" to "rectilinear" in the API to make this more clear?

I think it would be nice, if we could generate fields directly on pyvista meshes like you have proposed in #59 where these meshes are always recasted to unstructured meshes. Therefore it would also be nice to have a better handling of different generated fields in the class.

That'll work perfectly then - with a little bit a trickery, we cast the pyvista meshes to unstructured for GSTools to solve on them, then collect the result back on to structured grids

[LSchueler]

One more advantage of explicitly handling rectilinear grids, is the much easier handling of directional variograms. On all other grid types, you have to set an origin, a direction and at least one angle to span a cone, in which the directional variogram is estimated.
On rectilinear grids on the other hand, you can simply rotate your grid so that the axes align with the anisotropy without loss of generality and then estimate your directional variograms.

The structured representation of the fields saves some memory, which is lost again, if you would use a structure like numpy.meshgrid, which needs as much memory, as unstructured grids do.

What about the terms "unstructured" vs. "structured" or "rectilinear" for backwards compatibility?

[MuellerSeb]

@banesullivan GSTools now internally handles everthing as unstructured grids to simplify the code.
So casting points to an unstructured grid should be the default way of dealing with any structured meshes.

The interface for meshio meshes was also refactored. Maybe it's time to provide an interface in Field.mesh for pyvista meshes? 😉

[MuellerSeb]

@banesullivan I reopened #59 since we decided to make the mesh class a separate container that can be passed to Field.mesh. So I would be happy finishing this old PR.

[MuellerSeb]

Time to close this one.

• all meshes are handled as unstructured internally
• passing "rectilinear" is now interpreted as "structured" since the only keyword that is checked is "unstructured" (default)
• Field.mesh now accepts pyvista meshes