Understanding Meshes and GeoTables visualisation capabilities

General Usage
,
1

I’m currently testing Meshes.jl and GeoTables to see if they can be useful for my finite element analyses. Part of the logic still feels a bit weird to me (attaching data to a topology feels a bit restrictive, but that is probably due to my lack of understanding of the GIS background).

A quite common problem with finite element models is that you have elements with very different topologies that arise. Looking at the Meshes documentation, I figured everything would be working smoothly.

However, it seems that it is not possible to plot in a single pass point, surfacic and linear topologies. As a example I expected the code below to show triangles/quadrangles and segments.

using Meshes, GLMakie
points = rand(Meshes.Point3f,10)
connec = connect.([(1,8),(1,10),(1,2,6,5),(2,4,6),(4,3,5,6),(3,1,5)], Ngon)
mesh = SimpleMesh(points, connec)
viz(mesh,facetcolor=:red,showfacets=false, segmentsize=10)

but no segment seems to be plotted

image
image697×735 72.5 KB

A side note, if Ngon is not specified in the connect calls, the connectivity vector show Segments, Triangles and Quadrangles. If Ngon is specified, Segments are replaced by Ngon, but not other topologies.

Also, if pointsize is set, I would have expected some markers at the points of the mesh, but it seems that to plot some points one needs to call viz(vertices(mesh)), the pointsize works… but it does not seem to be able to customize the marker (e.g. 2d as in scatter or 3D as in meshscatter). In the picture below, the z-level looks off…
image

All in all, viz feels pretty stiff, so I was wondering if it is just me or if it is not intended for this kind of usage ?

Regarding GeoTables, if one has different topologies in a mesh, it does not seem to be possible to impose a specific color per topology type. Is that an intended behavior ?

Finally, it does not seem to be easy to manipulate coordinates contained in a mesh (and therefore in a geotable), I would like to display animated FE solutions and that looks quite complicated with Meshes, is it ? The main application here would be to support Observables but that does not seem to be working with viz

1 Like
2

Thanks for reading the documentation, it is not as extensive as we would like it, but it is a good resource. Consider checking the test suite to learn more.

Before I can address some of your questions, I need to correct some terminology. A topology is associated with a Mesh, not with the elements of the mesh. For example:

julia> grid = CartesianGrid(10, 10)
10×10 CartesianGrid{2,Float64}
  minimum: Point(0.0, 0.0)
  maximum: Point(10.0, 10.0)
  spacing: (1.0, 1.0)

julia> topology(grid)
10×10 GridTopology(aperiodic, aperiodic)

Thanks for raising this issue. We’ve been using another type of Domain to plot these heterogeneous collections of geometries, it is called the GeometrySet:

julia> seg = rand(Segment{3,Float64})
Segment{3,Float64}
├─ Point(0.3536386125565003, 0.6247989917632064, 0.4913770054195603)
└─ Point(0.2047399101445041, 0.2095146651312808, 0.6902567148054136)

julia> tri = rand(Triangle{3,Float64})
Triangle{3,Float64}
├─ Point(0.3843660557682562, 0.9357014403697765, 0.24859171505811706)
├─ Point(0.21751576012459806, 0.3317780757725377, 0.7491410074439363)
└─ Point(0.655459010564051, 0.02177834402229417, 0.5379742907171491)

julia> quad = rand(Quadrangle{3,Float64})
Quadrangle{3,Float64}
├─ Point(0.6712004755817398, 0.1435166211968616, 0.3061647251387838)
├─ Point(0.8984797891905305, 0.8362509698145396, 0.8126090669929493)
├─ Point(0.31728045916357084, 0.4167679007359413, 0.3229987385971217)
└─ Point(0.7674954065464198, 0.9407291962395742, 0.12801548362577664)

julia> gset = GeometrySet([seg, tri, quad])
3 GeometrySet{3,Float64}
├─ Segment((0.353639, 0.624799, 0.491377), (0.20474, 0.209515, 0.690257))
├─ Triangle((0.384366, 0.935701, 0.248592), (0.217516, 0.331778, 0.749141), (0.655459, 0.0217783, 0.537974))
└─ Quadrangle((0.6712, 0.143517, 0.306165), ..., (0.767495, 0.940729, 0.128015))

julia> viz(gset)

image

However, you found a bug in our Mesh visualization, which calls Makie.jl mesh, which in turn doesn’t handle meshes with geometries of mixed parametric dimension. We need to handle that properly in our preprocessing steps.

This is expected behavior. The docstring of connect explains the different methods. Ngon is assumed by default when the number of points is greater or equal to 3.

You are correct. You can either visualize the Mesh or its vertices. To visualize both, you need to use viz followed by viz!.

Can you elaborate a bit more on how we could make it more flexible? viz was designed to visualize Domain as explained in the GDSJL book. After you grasp the concept, you can visualize pretty much anything as a combination of multiple viz and viz! calls.

I am not sure I understood the questions, but the GeoTable exists to store data for each element of the Domain. The GDSJL book has various examples of viz calls with different colors, and we also provide a viewer for convenience.

What kinds of manipulations? We support various geometric transforms such as Translate, Rotate, Scale, Affine, … that you can apply directly to a Mesh, Domain or GeoTable. Consider watching our recent workshops to learn more:

We are already working on other coordinate transforms that are popular in GIS, but I think that is not your use case.

We have full support for Observables, maybe you are relying on an old version of the framework. Please share a MWE and we can take a look.

3

Thanks for your quick reply, I figured you would be the one to answer ^^.

Topology

Yes, it is one of the semantic differences we have between finite element models and GIS I guess.

GeometrySet

I saw you created an issue about conversion between SimpleMesh and GeometrySet. For FE computations it would be very practical to have a mesh object that supports different types of elements, because that is just the case all the time in that context.

Ngons in connect

That is what I understood from the docstring

Finally, the type PL can be ommitted. In this case, the indices are assumed to be connected as a Ngon or as a Segment.

as well but the code below does not seem to work accordingly, unless I’m misunderstanding something.

julia> connect.([(1,2),(1,2,3),(1,2,3,4)])
3-element Vector{Connectivity}:
 Segment(1, 2)
 Triangle(1, 2, 3)
 Quadrangle(1, 2, 3, 4)

julia> connect.([(1,2),(1,2,3),(1,2,3,4)],Ngon)
3-element Vector{Connectivity}:
 Ngon(1, 2)
 Triangle(1, 2, 3)
 Quadrangle(1, 2, 3, 4)

Viz flexibility

I think given the vast possibilities that seem to be offered by Meshes I was expecting to be able to separately configure all the different topologies in the mesh or set displayed, but again that may not be the purpose of viz.
Also, given that I would pretty difficult to make everybody happy given all the potential use cases, it would be nice to have a doc section about implementing additional viz-like methods on meshes and geotables. I know that looking at the code is useful but a piece of doc never hurts.

Geometry Transforms & Observables

I checked the transforms (by the way it looks like Affine is not documented on the website, unless I missed something). I the case of FE solutions, one often tries to plot something like deformed_point = original_point + scale * solution, something that does not seem to be supported at the moment. I looks similar to Affine but is different.

I have the latest version, it is more of a use case problem than a bug or something else.
The animation application would make scale an Observable. How you you use Observables with viz and a geotable ?
I tried with viz(<geotable>.geometry,color=geotable.<prop>) but that looks incompatible with observables.

4

I am not sure this is the case. Topology is a mathematical concept, which is mostly used in FEM and CAD applications. GIS rarely discusses topology.

That is already possible. The SimpleMesh type can store different geometry types. Can you please clarify why you think this is not the case?

In the second example you are trying to connect two points into a Ngon, but n-gons only exist with at least 3 vertices (area). With two points you can only create a Segment. We should problably throw an error in the Ngon constructor to require more than 2 points.

Most of our keyword options apply to individual geometries (color, pointsize, etc.). If some option is not accepting a vector of options, then please open an issue and we can see if it makes sense.

I am not understanding what you are missing. Can you please share a MWE with expected behavior?

Opened an issue: Add `Affine` to docs · Issue #771 · JuliaGeometry/Meshes.jl · GitHub

It is definetly supported. Also explained in the book.

We forward the observables to the Makie.jl builtin plots. Are you sure it is an issue in our viz recipe and not in the Makie.jl builtin plot function?

5

Well, in FEM, we also discuss element topology, not only mesh topology, but that is not a major concern as long as we all understand each other.

Well, in FEM, point, linear, planar or volumic elements are not the single element types that may be of interest. And you might want to add a specific representation of these other element types. A common use case would be e.g. RBEs (rigid body elements) these have a specific connectivity, so in a sense they could be converted to a group of segments I guess, but I was wondering if additional types of elements could be added easily.

I do not have an MWE right now, but as there are plenty of elements, or other features of meshes, like boundary conditions or else. These would require specialized methods and element-like types. Maybe in GIS you would like to plot drilling-well trajectories or else and have a dedicated type and viz for that. Is this kind of usage completely odd ?

Yes I think that would be better, because otherwise returning silently something that should not work is a bit surprising.

I went through chapter II of the book about Transforms. Even though plenty of transforms are described, I can’t figure how to build a pipeline for my use case. For instance, assuming I have a geotable built with.

gt = georef((X = <vector>, Y = <vector>, Z = <vector>), <mesh>)

Deforming the mesh with a given scaling factor <scale> would require a transform pipeline such as

gt |> Translate(scale*gt.X,sale*gt.Y,scale*gt.Z)

but that does not seem to be supported right now. Is there another transform for this kind of application that I missed ? Moreover, if one uses an Observable as scaling factor, I definitely can’t see how to buld a lifting pipeline that would take the initial geotable as input and allow for vizualization with animation. If you have an example of that, it would be perfect, but there is no reference to this in the book as far as I know.
Assuming that the lifting pipeline would return an observable gt with type Observable{GeoTable}, the problem for me is that if you want to plot the mesh (deformed or not) you need to call

viz(gt[].geometry)

which would defeat the purpose of the observable.

6

Can you share one example of such a type? Are you referring to high-order elements with vertices on faces and edges?

Again, I am not following the issue. It is hard to understand what you have in mind without a MWE. We already have “specialized methods and element-like types”. We already plot drillling-well trajectories, see CylindricalTrajectory: Trajectories · Meshes.jl

Opened an issue: Do not allow `Ngon` with N < 3 · Issue #776 · JuliaGeometry/Meshes.jl · GitHub

Can you please confirm that you don’t want the following?

pipe = Scale(s1, s2, s3) → Translate(t1, t2, t3)

In the pipeline above you first scale, and then translate the result by another set of constants. If you want to translate the points individually based on a vector field, that is not supported yet. It should be easy to cook something after we have a precise definition of the arguments you have in mind.

That is a great feature request. We didn’t try pipelines with Observable parameters yet. Adding to the TODO list: Observable transform pipelines · Issue #399 · JuliaEarth/GeoStats.jl · GitHub

7

Hum, not necessarily, but that could be one case indeed, although the Polyhedron type should support most of such element already. An example is the RBE I mention in my previous message, they generally link one point to a group of points.

To make it hopefully more understandable, let’s say that you want to locate say a sensor in your mesh. You may want to have a specific Seismometer element, that may have a vertex number as property (or a Point with coordinates), on which you could call a specialized viz method. For my current applications, I have a specialized Makie recipes for each feature of my meshes allowing to display multiple meaningful features, in different manners.
In addition, for some features, color information is not systematically the single or best choice. Arrows, contours or else might be of interest, then easing the implementation of other visualization strategies could be nice.

To go back to my first message, could you explain the choice of 2D (transparent?) markers for 3D vertices visualization ?

Well, that depends on what is in the si and ti. For that I can provide the following MWE

using Meshes, GeoTables
grid = CartesianGrid(10,10,10)
gt = georef((TX=rand(1000),TY=rand(1000),TZ=rand(1000)),grid)
gt |> Translate(scale*gt.TX,scale*gt.TY,scale*gt.TZ) #fails because not applied to geometry
gt.geometry |> Translate(scale*gt.TX,scale*gt.TY,scale*gt.TZ) #fails because of vectors
gt.geometry |> Translate(scale*gt.TX[1],scale*gt.TY[1],scale*gt.TZ[1]) #works

In the use case I describe, one scales some column information in the table, then adds that to the initial domain coordinates, so I do not think this is what you describe. I understand what you say as a scaling of the initial domain, then a translation by a scalar of the whole domain.

8

Can you elaborate on why you can’t reproduce the same behavior with view? If you have a GeoTable where this tag information is stored as a column for each element, you can easily geotable[geotable.tag .== 1, :] to get a SubGeoTable for visualization with viz. Does that cover your use case?

You have to basically think in terms of geospatial data science, where domain elements have information stored in the corresponding rows of the feature table. Here is a helper function:

function viztag(geotable, col, val)
  vals = geotable[:,col] 
  gt = geotable[vals .== val, :]
  viz(gt)
end

Alternatively, you can store the tags in an external vector and perform a view(mesh, inds) to get the corresponding SubDomain for visualization.

I am not sure I understand the question. If you do viz(vertices(mesh)) you don’t get any transparency by defaul.t If you add alpha = 0.5, then you get transparency. Is that the question?

That is not a valid MWE, it doesn’t have the definition of scale, nor it compiles because Translate doesn’t accept vectors as arguments.

9

I think there is a misunderstanding here. I understand that geotables allow to comb through the data, hence my interest. My question is : is it possible to (or more likely how hard is it) to add application-specific element types for storage in a GeoTable and for visualization with viz. I am not saying the current viz is not helpful, I am just saying it will never cover all potential visualizations of data over meshes. Therefore having a documented way of adding new elements, domains or viz methods would be nice. As an example, I just may want to plot arrows or pyramids based on a geotable with specific elements. Question is, how can one write the proper recipe for that, should I just add some Makie based recipes, or does it make sense to rely on viz or else…

See the following example

using Meshes, GLMakie
points = rand(Meshes.Point3f,10)
connec = connect.([(1,2,6,5),(2,4,6),(4,3,5,6),(3,1,5)], Ngon)
mesh = SimpleMesh(points, connec)
viz(vertices(mesh), pointsize=100)
meshscatter!([GLMakie.Point3(v.coords...) for v in vertices(mesh)])
scatter!([GLMakie.Point3(v.coords...) for v in vertices(mesh)], markersize=75)

This gives
image

FYI, my environment is completely up to date, but the points plotted with viz look transparent. I thought that could be just a z-level problem/artifact during the plotting, but the additional scatter does not suffer from that it seems. Regarding the markers themselves, I was wondering why scatter was preferred evend for 3D data plots, instead of meshscatter.

Yes, sorry, scale is just a scalar factor. I know this does not compile (hence the fail comments at the end of the line). I was just showing you what I tested, but did not work (because not supported) and and what worked but was not my use case.
To push this a bit further, I guess this could be summarized as : is it possible to define coordinate transformations on the domain based on information in the table of a geotable ? As far as I understood from the book, you propose coordinate transforms, that modify the domain with scalar transforms (1 for the whole mesh e.g. rotation or translation), but no coordinate transforms that can use the information in the table, do you ?

10

You can add a new geometry in Meshes.jl in a PR, and it will automatically work with the GeoTable, viz, etc.

You can already do that by creating Ray and Pyramid geometries for example. You can place all these geometries into a GeometrySet and that will be a valid Domain. The Mesh is a specific subtype of Domain that is described in terms of vertices, connectivities, etc.

julia> using Meshes

julia> r = rand(Ray{3,Float64})
Ray{3,Float64}
├─ p: Point(0.8544621114647082, 0.007290166206372128, 0.8953959324059896)
└─ v: Vec(0.9568652328665942, 0.611726781018889, 0.27683459777111596)

julia> p = rand(Pyramid{3,Float64})
Pyramid{3,Float64}
├─ Point(0.5245605139763986, 0.31180747839374345, 0.8006794534270958)
├─ Point(0.17999965079093627, 0.2960387752699939, 0.23121858402483375)
├─ Point(0.16164742029403556, 0.5896417963031546, 0.0005634795956463989)
├─ Point(0.5200829045516021, 0.23759556572538754, 0.36635011529535155)
└─ Point(0.8140821418410266, 0.522107532502354, 0.5031798558753174)

julia> d = GeometrySet([r, p])
2 GeometrySet{3,Float64}
├─ Ray(p: (0.854462, 0.00729017, 0.895396), v: (0.956865, 0.611727, 0.276835))
└─ Pyramid((0.524561, 0.311807, 0.800679), ..., (0.814082, 0.522108, 0.50318))

julia> import GLMakie as Mke

julia> viz(rand(Ray{3,Float64}, 100))

image

Here is the full recipe to plot a PointSet, which is constructed when you ask to viz the vertices of the mesh:

We simply call Makie.scatter with the overdraw=true keyword to address issues with large point sets, see Hide stroke in 3D scatter using Makie.jl.

The CoordinateTransform is a special case of GeometricTransform that simply applies the same function to each point of the domain individually. These transforms are defined in terms of parameters, they do not live in the columns of the GeoTable. If you want to design a transform that takes both the values and the domain of the GeoTable, then you are after a general TableTransform. See the hierarchy of transforms here: Transforms · GeoStats.jl

I think the easiest way forward is to start adding missing geometries in Meshes.jl, then trying to visualize them. After that, we can come up with higher-level transforms that meet your requirements, and a find a good home for them. Meshes.jl only defines GeometricTransform that act on the domain only, they don’t have access to the table of features.

1 Like
11

I have one more question @juliohm, I tried to find where the scene settings were defined, but I did not find them. The lightning options in viz lead to weird results with the mesh I’m testing with, and I would like to revert to the Makie defaults.

Even on a simple cube it looks slightly funny to my taste :slight_smile:

using Meshes, GLMakie
f =  Figure()
# Meshes cubic grid
l = LScene(f[1,1])
grid = CartesianGrid(10,10,10)
viz!(grid)
# GLMakie cubic grid
l = LScene(f[1,2])

rectmesh = Rect3(GLMakie.Point3f(-0.5), GLMakie.Vec3f(1))
mesh!(rectmesh)

image

12

I think you are plotting different things. On the left you have a CartesianGrid, on the right you have a Rect. Try to viz a Meshes.Box or Meshes.Hexahedron, you can also enable the option showfacets=true to highlight the edges.

If some additional option is required, we are happy to forward it to the underlying Makie function call. The viz is nothing more than a recipe that does preprocessing and calls Makie. It saves you from a lot of low-level stuff.

13

Here is the link for the recipe you called, a 3D Cartesian grid:

If you prefer to discuss this further in a quick chat, we are on Zulip. There it is easier to iterate the details:

https://julialang.zulipchat.com/#narrow/stream/275558-meshes.2Ejl

14

You can already apply transform pipelines with Mke.Observable parameters. This is a MWE with real-time translation of grids:

using GeoStats # load stack

import GLMakie as Mke # load Makie backend

grid = CartesianGrid(10, 10) # domain of interest

vec = Mke.Observable((2, 2)) # translation vector

tgrid = Mke.@lift grid |> Translate($vec) # apply translation

viz(tgrid, showsegments=true) # visualize results

vec[] = (-2, -2) # update translation vector