Can I put non-primitive geometries into a GeometrySet?

General Usage
1

If I put a SimpleMesh geometry into a GeometrySet, I get a StackOverflowError:

using Meshes

# test using primitive geometries
s1 = Meshes.Sphere((0, 0, 0), 1)
s2 = Meshes.Sphere((1, 1, 1), 1)
gs1 = GeometrySet([s1, s2]) # works fine

# test using Meshes
s3 = discretize(s1, RegularDiscretization(10, 10))
s4 = discretize(s2, RegularDiscretization(10, 10))
gs2 = GeometrySet([s3, s4]) # throws StackOverflowError

StackOverflowError:

Stacktrace:
  [1] _stable_typeof
    @ ./operators.jl:910 [inlined]
  [2] Base.Fix2(f::typeof(isequal), x::Type)
    @ Base ./operators.jl:1121
  [3] isequal
    @ ./operators.jl:1136 [inlined]
  [4] allequal
    @ ./set.jl:533 [inlined]
  [5] GeometrySet(geoms::Vector{SimpleMesh{๐”ผ{3}, Cartesian3D{NoDatum, Quantity{Float64, ๐‹, Unitful.FreeUnits{(m,), ๐‹, nothing}}}, Vector{Meshes.Point{๐”ผ{3}, Cartesian3D{NoDatum, Quantity{Float64, ๐‹, Unitful.FreeUnits{(m,), ๐‹, nothing}}}}}, SimpleTopology{Connectivity}}})
    @ Meshes ~/.julia/packages/Meshes/fFipW/src/domains/sets.jl:29
  [6] GeometrySet(geoms::Vector{SimpleMesh{๐”ผ{3}, Cartesian3D{NoDatum, Quantity{Float64, ๐‹, Unitful.FreeUnits{(m,), ๐‹, nothing}}}, Vector{Meshes.Point{๐”ผ{3}, Cartesian3D{NoDatum, Quantity{Float64, ๐‹, Unitful.FreeUnits{(m,), ๐‹, nothing}}}}}, SimpleTopology{Connectivity}}}) (repeats 74533 times)
    @ Meshes ~/.julia/packages/Meshes/fFipW/src/domains/sets.jl:35
  [7] eval
    @ ./boot.jl:370 [inlined]
  [8] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1903
...

Thoughts, suggestions? Iโ€™m new to Meshes.jl, so very possibly doing something stupid. I am headed in this direction (meshes in GeometrySets) because ultimately I want to save multiple geometries into an OBJ or PLY file and it seems like using a geotable and GeoIO would be a convenient way to do that.

Many thanks!

Iโ€™m using Julia 1.9.2 and Meshes.jl 0.51.19

2

In Meshes.jl we differentiate between Geometry and Domain. The Domain is a collection of geometries. The GeometrySet is a Domain where the geometries know nothing about each other.

Can I put non-primitive geometries into a GeometrySet?

You can store any Geometry in a GeometrySet, including Primitive and Polytope geometries.

Did you have a chance to read the GDSJL book?

https://juliaearth.github.io/geospatial-data-science-with-julia

3

The GeoTable stores a single Domain in the :geometry column. GeoIO.jl can only load/save one GeoTable at a time. If you need to save two meshes that are the result of discretization of geometries, you can try to merge them into a single mesh.

Based on your example:

mesh1 = discretize(s1, RegularDiscretization(10, 10))
mesh2 = discretize(s2, RegularDiscretization(10, 10))

mesh = merge(mesh1, mesh2)

geotable = georef(nothing, mesh) # without attributes

GeoIO.save("file.ply", geotable)
4

Thank you very much, Julio. Very helpful and timely replies. I was indeed being stupid and was loose in my understanding of Domain (even though I had read GDSJLโ€ฆ too quickly).

I am unfortunately having some trouble with the practicalities of merging. If I create meshes from a box and a sphere and merge them, the new mesh seems to have all the vertices of the combination, but only the connectivity from the box. Hereโ€™s an example:

box1_mesh = discretize(Meshes.Box((1, 1, 1), (2, 2, 2)), RegularDiscretization(1, 1))
sphere1_mesh = discretize(Meshes.Sphere((2, 2, 2), 0.5), RegularDiscretization(4, 4))
# one box and one sphere do not merge as expected (by me)
oneeach = merge(box1_mesh, sphere1_mesh)
1 SimpleMesh
  30 vertices
  โ”œโ”€ Point(x: 1.0 m, y: 1.0 m, z: 1.0 m)
  โ”œโ”€ Point(x: 2.0 m, y: 1.0 m, z: 1.0 m)
  โ”œโ”€ Point(x: 1.0 m, y: 2.0 m, z: 1.0 m)
  โ”œโ”€ Point(x: 2.0 m, y: 2.0 m, z: 1.0 m)
  โ”œโ”€ Point(x: 1.0 m, y: 1.0 m, z: 2.0 m)
  โ‹ฎ
  โ”œโ”€ Point(x: 2.0 m, y: 1.5 m, z: 2.0 m)
  โ”œโ”€ Point(x: 2.0 m, y: 1.5669872981077806 m, z: 1.75 m)
  โ”œโ”€ Point(x: 2.0 m, y: 1.75 m, z: 1.5669872981077806 m)
  โ”œโ”€ Point(x: 2.0 m, y: 2.0 m, z: 2.5 m)
  โ””โ”€ Point(x: 2.0 m, y: 2.0 m, z: 1.5 m)
  1 elements
  โ””โ”€ Hexahedron(1, 2, 4, 3, 5, 6, 8, 7)

Merging two box meshes or two sphere meshes works as I expected. What am I missing this time? :smile: Thanks for your time!

5

Hi @Jfschue , that was a tricky one. Notice that the first mesh is made of 3D geometries and the second is made of 2D geometries. When you merge the two, the SimpleMesh constructor will take all these connectivities into account, but will only interpret the geometries with highest dimension as the elements.

You can see that all connectivities are stored in the underlying topology:

topology(oneeach)

SimpleTopology{Connectivity}(Connectivity[Hexahedron(1, 2, 4, 3, 5, 6, 8, 7), Quadrangle(9, 10, 15, 14), Quadrangle(10, 11, 16, 15), Quadrangle(11, 12, 17, 16), Quadrangle(12, 13, 18, 17), Quadrangle(14, 15, 20, 19), Quadrangle(15, 16, 21, 20), Quadrangle(16, 17, 22, 21), Quadrangle(17, 18, 23, 22), Quadrangle(19, 20, 25, 24)  โ€ฆ  Quadrangle(26, 27, 12, 11), Quadrangle(27, 28, 13, 12), Triangle(29, 9, 14), Triangle(29, 14, 19), Triangle(29, 19, 24), Triangle(29, 24, 9), Triangle(30, 18, 13), Triangle(30, 23, 18), Triangle(30, 28, 23), Triangle(30, 13, 28)], [3, 2, 2, 2, 2, 2, 2, 2, 2, 2  โ€ฆ  2, 2, 2, 2, 2, 2, 2, 2, 2, 2], [1])

It is saying that you have 1 Hexahedron (the single element) and a couple of Quadrangle and Triangle of lower parametric dimension.

The GeoIO.jl module currently loads/saves element and vertex data. It ignores geometries with intermediate dimension.

6

Unbounded recursion is a bug, independently of what specifically caused it. Ideally the constructor would be coded in such a manner to make infinite recursion impossible (by throwing explicitly or, better, relying on dispatch to throw a MethodError). So, in principle, this is a package bug, to be pedantic.

1 Like
7

And we can find many other bugs that way when a constructor is called with unexpected arguments. PRs are welcome to try to improve the error message and avoid the stack overflow. Not in our priority list to fix bugs caused by incorrect usage.

1 Like
8

Thanks for the explanation. Good to know whatโ€™s going on. This detail is not a limitation for my actual use case, but Iโ€™m glad it came up here. I learned something new.

1 Like
9

You make an excellent point @nsajko. There are always opportunities to inform [new] users when things go unexpectedly. I read the following somewhereโ€ฆ

Geometric processing doesnโ€™t need to be complicated. It should be fun and read like math. If it feels โ€œcomputer sciencyโ€, that is a limitation of the software and programming language.

Thanks for the contribution!

1 Like