Issubset gives error for this Hexahedron

[raman-maker]

In this code m ⊆ m2 should output either true or false. See this https://discourse.julialang.org/t/issubset-gives-error-for-this-hexahedron/133346

using Meshes

m = Hexahedron((0.1547376288700546, 0.05921552849426126, 1.7379087983473205), (0.17193069874450512, 0.05921552849426126, 1.7379087983473205), (0.17193069874450512, 0.06579503166029028, 1.7379087983473205), (0.1547376288700546, 0.06579503166029028, 1.7379087983473205),(0.1547376288700546, 0.05921552849426126, 1.754088044166565),(0.17193069874450512, 0.05921552849426126, 1.754088044166565), (0.17193069874450512, 0.06579503166029028, 1.754088044166565), (0.1547376288700546, 0.06579503166029028, 1.754088044166565))
m2 = Hexahedron((0.0, 0.0, 1.600000023841858), (0.0, 0.0, 1.754088044166565),(0.1719306987445051, 0.0, 1.7456416303231632),  (0.15682743121429854, 0.0, 1.5922955859741212), (0.0, 0.0, 1.600000023841858), (0.0, 0.0, 1.754088044166565), (0.15884325286141418, 0.06579503166029028, 1.7456416303231632), (0.1448896531793745, 0.06001526124941982, 1.5922955859741212))

m ⊆ m2

But it gives error.

ERROR: MethodError: no method matching ntuple(::Meshes.var"#apply##14#apply##15"{Meshes.Point{…}}, ::Quantity{Float64, 𝐋, Unitful.FreeUnits{…}})
The function ntuple exists, but no method is defined for this combination of argument types.

The problem arises in issubset.jl.

[JoshuaLampert]

For reference, the whole stacktrace is

ERROR: MethodError: no method matching ntuple(::Meshes.var"#apply##14#apply##15"{Point{…}}, ::Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{…}})
The function `ntuple` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  ntuple(::Any, ::Val{3})
   @ Base ntuple.jl:52
  ntuple(::Any, ::Val{2})
   @ Base ntuple.jl:51
  ntuple(::Any, ::Val{1})
   @ Base ntuple.jl:50
  ...

Stacktrace:
  [1] apply(::Repair{8}, poly::Quadrangle{𝔼{3}, CoordRefSystems.Cartesian3D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}})
    @ Meshes ~/.julia/dev/Meshes/src/transforms/repair.jl:124
  [2] Transform
    @ ~/.julia/packages/TransformsBase/w6gll/src/interface.jl:124 [inlined]
  [3] |>
    @ ./operators.jl:972 [inlined]
  [4] discretize(polygon::Quadrangle{𝔼{3}, CoordRefSystems.Cartesian3D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}, method::DehnTriangulation)
    @ Meshes ~/.julia/dev/Meshes/src/discretization.jl:104
  [5] simplexify
    @ ~/.julia/dev/Meshes/src/discretization.jl:208 [inlined]
  [6] _appendinds!(ginds::Vector{Vector{Int64}}, connec::Connectivity{Quadrangle, 4}, points::Vector{Point{𝔼{3}, CoordRefSystems.Cartesian3D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}})
    @ Meshes ~/.julia/dev/Meshes/src/discretization.jl:259
  [7] _simplexify(points::Vector{Point{𝔼{3}, CoordRefSystems.Cartesian3D{CoordRefSystems.NoDatum, Unitful.Quantity{…}}}}, connec::Base.Generator{UnitRange{Int64}, Meshes.var"#elements##0#elements##1"{SimpleTopology{Connectivity{…}}}})
    @ Meshes ~/.julia/dev/Meshes/src/discretization.jl:237
  [8] simplexify
    @ ~/.julia/dev/Meshes/src/discretization.jl:223 [inlined]
  [9] winding(points::Vector{Point{𝔼{…}, CoordRefSystems.Cartesian3D{…}}}, mesh::SimpleMesh{𝔼{3}, CoordRefSystems.Cartesian3D{CoordRefSystems.NoDatum, Unitful.Quantity{…}}, Vector{Point{…}}, SimpleTopology{Connectivity{…}}})
    @ Meshes ~/.julia/dev/Meshes/src/winding.jl:46
 [10] sidewithinbox
    @ ~/.julia/dev/Meshes/src/sideof.jl:201 [inlined]
 [11] sideof(points::Tuple{Point{𝔼{…}, CoordRefSystems.Cartesian3D{…}}}, object::SimpleMesh{𝔼{3}, CoordRefSystems.Cartesian3D{CoordRefSystems.NoDatum, Unitful.Quantity{…}}, Vector{Point{…}}, SimpleTopology{Connectivity{…}}})
    @ Meshes ~/.julia/dev/Meshes/src/sideof.jl:195
 [12] sideof
    @ ~/.julia/dev/Meshes/src/sideof.jl:182 [inlined]
 [13] in
    @ ~/.julia/dev/Meshes/src/predicates/in.jl:10 [inlined]
 [14] Fix
    @ ./operators.jl:1194 [inlined]
 [15] #246
    @ ~/.julia/packages/StaticArrays/DsPgf/src/mapreduce.jl:305 [inlined]
 [16] macro expansion
    @ ~/.julia/packages/StaticArrays/DsPgf/src/mapreduce.jl:180 [inlined]
 [17] _mapfoldl(f::StaticArrays.var"#246#247"{Base.Fix2{typeof(in), Hexahedron{…}}}, op::typeof(&), dims::Colon, init::Bool, ::StaticArraysCore.Size{(8,)}, a::StaticArraysCore.SVector{8, Point{𝔼{…}, CoordRefSystems.Cartesian3D{…}}})
    @ StaticArrays ~/.julia/packages/StaticArrays/DsPgf/src/mapreduce.jl:155
 [18] _mapreduce(::Function, ::Function, ::Function, ::Bool, ::StaticArraysCore.Size{(8,)}, ::StaticArraysCore.SVector{8, Point{𝔼{3}, CoordRefSystems.Cartesian3D{CoordRefSystems.NoDatum, Unitful.Quantity{…}}}})
    @ StaticArrays ~/.julia/packages/StaticArrays/DsPgf/src/mapreduce.jl:153
 [19] all
    @ ~/.julia/packages/StaticArrays/DsPgf/src/mapreduce.jl:305 [inlined]
 [20] issubset
    @ ~/.julia/dev/Meshes/src/predicates/issubset.jl:38 [inlined]
 [21] issubset(g₁::Hexahedron{CoordRefSystems.Cartesian3D{CoordRefSystems.NoDatum, Unitful.Quantity{…}}, 𝔼{3}}, g₂::Hexahedron{CoordRefSystems.Cartesian3D{CoordRefSystems.NoDatum, Unitful.Quantity{…}}, 𝔼{3}})
    @ Meshes ~/.julia/dev/Meshes/src/predicates/issubset.jl:14
 [22] top-level scope
    @ REPL[8]:1
Some type information was truncated. Use `show(err)` to see complete types.

Reducing this even further leads to

julia> elem = Quadrangle(Point(0.0, 0.0, 1.754088044166565), Point(0.0, 0.0, 1.600000023841858),
                         Point(0.0, 0.0, 1.600000023841858), Point(0.0, 0.0, 1.754088044166565))
Quadrangle
├─ Point(x: 0.0 m, y: 0.0 m, z: 1.754088044166565 m)
├─ Point(x: 0.0 m, y: 0.0 m, z: 1.600000023841858 m)
├─ Point(x: 0.0 m, y: 0.0 m, z: 1.600000023841858 m)
└─ Point(x: 0.0 m, y: 0.0 m, z: 1.754088044166565 m)

julia> simplexify(elem)
ERROR: MethodError: no method matching ntuple(::Meshes.var"#apply##14#apply##15"{Point{…}}, ::Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{…}})
The function `ntuple` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  ntuple(::Any, ::Val{3})
   @ Base ntuple.jl:52
  ntuple(::Any, ::Val{2})
   @ Base ntuple.jl:51
  ntuple(::Any, ::Val{1})
   @ Base ntuple.jl:50
  ...

Stacktrace:
 [1] apply(::Repair{8}, poly::Quadrangle{𝔼{3}, CoordRefSystems.Cartesian3D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}})
   @ Meshes ~/.julia/dev/Meshes/src/transforms/repair.jl:124
 [2] Transform
   @ ~/.julia/packages/TransformsBase/w6gll/src/interface.jl:124 [inlined]
 [3] |>
   @ ./operators.jl:972 [inlined]
 [4] discretize(polygon::Quadrangle{𝔼{3}, CoordRefSystems.Cartesian3D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}, method::DehnTriangulation)
   @ Meshes ~/.julia/dev/Meshes/src/discretization.jl:104
 [5] simplexify(poly::Quadrangle{𝔼{3}, CoordRefSystems.Cartesian3D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}})
   @ Meshes ~/.julia/dev/Meshes/src/discretization.jl:208
 [6] top-level scope
   @ REPL[4]:1
Some type information was truncated. Use `show(err)` to see complete types.

Here, elem is one element in boundary(m2). So it looks like this Quadrangle is degenerated, which causes this problem.

[raman-maker]

Here, elem is one element in boundary(m2). So it looks like this Quadrangle is degenerated, which causes this problem.

I think my data is correct. I tried code with other files and i got same error. So error is in Meshes.jl implementation. see this new data of hexahedrons.

using Meshes

h = Hexahedron((0.05054647164016874, 0.009861127699934922, 1.140005387173362), (0.056162746266854156, 0.009861127699934922, 1.140005387173362), (0.056162746266854156, 0.010956808555483246, 1.140005387173362), (0.05054647164016874, 0.010956808555483246, 1.140005387173362),  (0.05054647164016874, 0.009861127699934922, 1.1445976495742798),(0.056162746266854156, 0.009861127699934922, 1.1445976495742798), (0.056162746266854156, 0.010956808555483246, 1.1445976495742798), (0.05054647164016874, 0.010956808555483246, 1.1445976495742798))
h2 = Hexahedron((0.0, 0.0, 1.100000023841858),  (0.0, 0.0, 1.1445976495742798), (0.056162746266854156, 0.0, 1.1432189315011936), (0.05397444443079348, 0.0, 1.098675025565101), (0.0, 0.0, 1.100000023841858), (0.0, 0.0, 1.1445976495742798), (0.05508359478568487, 0.010956808555483246, 1.1432189315011936), (0.05293734055812983, 0.010529892033890714, 1.098675025565101))

h ⊆ h2

What should i do?

[juliohm]

Thank you @JoshuaLampert for looking into it. I am always skeptical about the issues and Discourse threads opened by this user.

It doesn't matter how many times we try to explain the basics. I am starting to feel the account is fake or intentionally created to expose "bugs" that don't exist.

[raman-maker]

Hi @juliohm This data is from simulation and is not fake. I am human and not robot. If you like i can share data with you. And nobody pays me to intentionally create bugs.

[JoshuaLampert]

To elaborate a bit more on this. Your m2 above (and also your h2) are no Hexahedrons (or call them degenerate Hexahedrons), but a triangular prism or wedge. If you look at a plot

you'll see that it only has 5 faces instead of 6 due to how you pick the points (there are duplicates). This leads to the degenerate Quadrangles in the mesh describing the boundary:

julia> mesh = boundary(m2)
6 SimpleMesh
  8 vertices
  ├─ Point(x: 0.0 m, y: 0.0 m, z: 1.600000023841858 m)
  ├─ Point(x: 0.0 m, y: 0.0 m, z: 1.754088044166565 m)
  ├─ Point(x: 0.1719306987445051 m, y: 0.0 m, z: 1.7456416303231632 m)
  ├─ Point(x: 0.15682743121429854 m, y: 0.0 m, z: 1.5922955859741212 m)
  ├─ Point(x: 0.0 m, y: 0.0 m, z: 1.600000023841858 m)
  ├─ Point(x: 0.0 m, y: 0.0 m, z: 1.754088044166565 m)
  ├─ Point(x: 0.15884325286141418 m, y: 0.06579503166029028 m, z: 1.7456416303231632 m)
  └─ Point(x: 0.1448896531793745 m, y: 0.06001526124941982 m, z: 1.5922955859741212 m)
  6 elements
  ├─ Quadrangle(4, 3, 2, 1)
  ├─ Quadrangle(6, 5, 1, 2)
  ├─ Quadrangle(3, 7, 6, 2)
  ├─ Quadrangle(4, 8, 7, 3)
  ├─ Quadrangle(1, 5, 8, 4)
  └─ Quadrangle(6, 7, 8, 5)

julia> points = vertices(mesh)
8-element Vector{Meshes.Point{𝔼{3}, CoordRefSystems.Cartesian3D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}}:
 Point(x: 0.0 m, y: 0.0 m, z: 1.600000023841858 m)
 Point(x: 0.0 m, y: 0.0 m, z: 1.754088044166565 m)
 Point(x: 0.1719306987445051 m, y: 0.0 m, z: 1.7456416303231632 m)
 Point(x: 0.15682743121429854 m, y: 0.0 m, z: 1.5922955859741212 m)
 Point(x: 0.0 m, y: 0.0 m, z: 1.600000023841858 m)
 Point(x: 0.0 m, y: 0.0 m, z: 1.754088044166565 m)
 Point(x: 0.15884325286141418 m, y: 0.06579503166029028 m, z: 1.7456416303231632 m)
 Point(x: 0.1448896531793745 m, y: 0.06001526124941982 m, z: 1.5922955859741212 m)

julia> connec = elements(topology(mesh))
Base.Generator{UnitRange{Int64}, Meshes.var"#elements##0#elements##1"{SimpleTopology{Connectivity{Quadrangle, 4}}}}(Meshes.var"#elements##0#elements##1"{SimpleTopology{Connectivity{Quadrangle, 4}}}(SimpleTopology{Connectivity{Quadrangle, 4}}(Connectivity{Quadrangle, 4}[Quadrangle(4, 3, 2, 1), Quadrangle(6, 5, 1, 2), Quadrangle(3, 7, 6, 2), Quadrangle(4, 8, 7, 3), Quadrangle(1, 5, 8, 4), Quadrangle(6, 7, 8, 5)], [2, 2, 2, 2, 2, 2], [1, 2, 3, 4, 5, 6])), 1:6)

julia> c = collect(connec)[2]
Quadrangle(6, 5, 1, 2)

julia> elem = materialize(c, points)
Quadrangle
├─ Point(x: 0.0 m, y: 0.0 m, z: 1.754088044166565 m)
├─ Point(x: 0.0 m, y: 0.0 m, z: 1.600000023841858 m)
├─ Point(x: 0.0 m, y: 0.0 m, z: 1.600000023841858 m)
└─ Point(x: 0.0 m, y: 0.0 m, z: 1.754088044166565 m)

You could try to use a Wedge instead by removing the doubled vertices from above, e.g.,

julia> m2 = Wedge((0.0, 0.0, 1.754088044166565),(0.1719306987445051, 0.0, 1.7456416303231632),  (0.15682743121429854, 0.0, 1.5922955859741212), (0.0, 0.0, 1.600000023841858), (0.15884325286141418, 0.06579503166029028, 1.7456416303231632), (0.1448896531793745, 0.06001526124941982, 1.5922955859741212))
Wedge
├─ Point(x: 0.0 m, y: 0.0 m, z: 1.754088044166565 m)
├─ Point(x: 0.1719306987445051 m, y: 0.0 m, z: 1.7456416303231632 m)
├─ Point(x: 0.15682743121429854 m, y: 0.0 m, z: 1.5922955859741212 m)
├─ Point(x: 0.0 m, y: 0.0 m, z: 1.600000023841858 m)
├─ Point(x: 0.15884325286141418 m, y: 0.06579503166029028 m, z: 1.7456416303231632 m)
└─ Point(x: 0.1448896531793745 m, y: 0.06001526124941982 m, z: 1.5922955859741212 m)

If you vizualize this Wedge, you'll see it's the same as the degenerate Hexahedron from above. However, Meshes.jl currently does not implement the function isconvex for a Wedge (similar to #1205), which means you will still not be able to call m ⊆ m2.

[raman-maker]

This gives true. I think Hexahedron struct should be updated may be using inner constructor that checks redundancy of vertices.

using Meshes
m2 = Hexahedron((0.0, 0.0, 1.600000023841858), (0.0, 0.0, 1.754088044166565),(0.1719306987445051, 0.0, 1.7456416303231632),  (0.15682743121429854, 0.0, 1.5922955859741212), (0.0, 0.0, 1.600000023841858), (0.0, 0.0, 1.754088044166565), (0.15884325286141418, 0.06579503166029028, 1.7456416303231632), (0.1448896531793745, 0.06001526124941982, 1.5922955859741212))
m2 isa Hexahedron

Should i open issue for isconvex for a Wedge? Will opening new issue for isconvex for a Wedge hurt @juliohm ?

[JoshuaLampert]

Of course this will give true because m2 is of type Hexahedron. What I meant with m2 is not a hexahedron was that geometrically it is not a hexahedron. What we could do in Meshes.jl is to check on creation of a Hexahedron if there are duplicated vertices and immediately error in this case, but I doubt this is really helpful and would be something we want.

[KristofferC]

Even for degenerate shapes, a better error message than

ERROR: MethodError: no method matching ntuple(

is surely desirable, or?

[JoshuaLampert]

It's worth considering to improve error messages. However, I'm not sure where the improved error message should live. Do we want to already disallow creating degenerate geometries or should we rather throw some descriptive errors when one applies an algorithm on a degenerate geometry, which assumes non-degeneracy. If we want this, I feel like this needs to be more thought-out.

[juliohm]

We can certainly consider better error messages when users do their homework to use the package as expected. It is an endless effort otherwise to intercept all possible inputs from users to throw errors and slow things down. In this case I believe the error message should be thrown by the specific function that assumes non-degenerate hexahedra. Alternatively, someone who really needs the degenerative case can jump in and contribute a generalization of the algorithm. Em qua., 22 de out. de 2025, 12:22, Joshua Lampert ***@***.***> escreveu:

…

*JoshuaLampert* left a comment (JuliaGeometry/Meshes.jl#1260) <#1260 (comment)> It's worth considering to improve error messages. However, I'm not sure where the improved error message should live. Do we want to already disallow creating degenerate geometries or should we rather throw some descriptive errors when one applies an algorithm on a degenerate geometry, which assumes non-degeneracy. If we want this, I feel like this needs to be more thought-out. — Reply to this email directly, view it on GitHub <#1260 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAZQW3K74SSSLRBIXYCIWCL3Y6OLNAVCNFSM6AAAAACJ4MCGZKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZTIMZSHE4DMMJYGE> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[raman-maker]

• My problem was this long code. Spherical grid hexahedron el don't have duplicate points but sel have duplicate points. Duplication occurs when i convert them to Cartesian coordinates.
• I want to find a bounding box for el or sel. I will then create a RegularGrid on this bounding box and keep only those RegulGrid elements of box that are completely inside el or sel.
• Since these RegularGrid elements are cuboids so i can get their 3D Makie.volume!() plot and stitch together to get final volume plot that is not possible for irregular hexahedron or prisms..

using Meshes, CoordRefSystems ,GLMakie, LinearAlgebra, Unitful
 
r = Float32[1.1, 1.1445976, 1.1910034, 1.2392906, 1.2895356, 1.3418176, 1.3962194, 1.4528267, 1.5117291, 1.5730196, 1.636795, 1.7031561, 1.7722077, 1.8440589, 1.9188232, 1.9966186, 2.077568, 2.1617997, 2.2494462, 2.3406463, 2.4355438, 2.5342886, 2.6370373]
θ = Float32[0.0, 0.049087387, 0.09817477, 0.14726216, 0.19634955]
ϕ = Float32[0.0, 0.19634955, 0.3926991, 0.5890486, 0.7853982, 0.9817477, 1.1780972, 1.3744467, 1.5707964, 1.7671459, 1.9634954, 2.1598449, 2.3561945, 2.552544, 2.7488935, 2.9452431, 3.1415927]

g = RectilinearGrid{𝔼,typeof(Spherical(0,0,0))}(r, θ, ϕ)
el = element(g,1)
vc = coords.(vertices(el))
cvc = convert.(Cartesian, vc)
points = [Meshes.Point(Float64.(ustrip.((c.x, c.y, c.z)))...) for c in cvc]
sel =Hexahedron(tuple(points...))
bb = Meshes.boundingbox(sel)

# Collect voxels (cuboids) inside the hexahedron
dbb = RegularGrid(bb.min, bb.max, dims = (10, 10, 10))
cuboids = []
for cl in nelements(dbb)
	c = element(dbb,cl)
	c ⊆ sel
	if c ⊆ sel     # keep cuboid if inside element
		push!(cuboids, cl)
	end
end

julia> el
Hexahedron
├─ Point(r: 1.100000023841858 m, θ: 0.0 rad, ϕ: 0.0 rad)
├─ Point(r: 1.1445976495742798 m, θ: 0.0 rad, ϕ: 0.0 rad)
├─ Point(r: 1.1445976495742798 m, θ: 0.04908738657832146 rad, ϕ: 0.0 rad)
├─ Point(r: 1.100000023841858 m, θ: 0.04908738657832146 rad, ϕ: 0.0 rad)
├─ Point(r: 1.100000023841858 m, θ: 0.0 rad, ϕ: 0.19634954631328583 rad)
├─ Point(r: 1.1445976495742798 m, θ: 0.0 rad, ϕ: 0.19634954631328583 rad)
├─ Point(r: 1.1445976495742798 m, θ: 0.04908738657832146 rad, ϕ: 0.19634954631328583 rad)
└─ Point(r: 1.100000023841858 m, θ: 0.04908738657832146 rad, ϕ: 0.19634954631328583 rad)

julia> sel
Hexahedron
├─ Point(x: 0.0 m, y: 0.0 m, z: 1.100000023841858 m)
├─ Point(x: 0.0 m, y: 0.0 m, z: 1.1445976495742798 m)
├─ Point(x: 0.056162746266854156 m, y: 0.0 m, z: 1.1432189315011936 m)
├─ Point(x: 0.05397444443079348 m, y: 0.0 m, z: 1.098675025565101 m)
├─ Point(x: 0.0 m, y: 0.0 m, z: 1.100000023841858 m)
├─ Point(x: 0.0 m, y: 0.0 m, z: 1.1445976495742798 m)
├─ Point(x: 0.05508359478568487 m, y: 0.010956808555483246 m, z: 1.1432189315011936 m)
└─ Point(x: 0.05293734055812983 m, y: 0.010529892033890714 m, z: 1.098675025565101 m)

[juliohm]

And that is why I say this user is the professional spammer. He copied/pasted this same example dozens of times on different Discourse threads, Makie.jl issues and now here. Please ignore it. Em qua., 22 de out. de 2025, 14:34, Raman Kumar ***@***.***> escreveu:

…

*raman-maker* left a comment (JuliaGeometry/Meshes.jl#1260) <#1260 (comment)> My problem was this long code. Spherical grid hexahedron el don't have duplicate points. Duplication occurs when i convert them to Cartesian coordinates. using Meshes, CoordRefSystems ,GLMakie, LinearAlgebra, Unitful r = Float32[1.1, 1.1445976, 1.1910034, 1.2392906, 1.2895356, 1.3418176, 1.3962194, 1.4528267, 1.5117291, 1.5730196, 1.636795, 1.7031561, 1.7722077, 1.8440589, 1.9188232, 1.9966186, 2.077568, 2.1617997, 2.2494462, 2.3406463, 2.4355438, 2.5342886, 2.6370373] θ = Float32[0.0, 0.049087387, 0.09817477, 0.14726216, 0.19634955] ϕ = Float32[0.0, 0.19634955, 0.3926991, 0.5890486, 0.7853982, 0.9817477, 1.1780972, 1.3744467, 1.5707964, 1.7671459, 1.9634954, 2.1598449, 2.3561945, 2.552544, 2.7488935, 2.9452431, 3.1415927] g = RectilinearGrid{𝔼,typeof(Spherical(0,0,0))}(r, θ, ϕ) el = element(g,1) vc = coords.(vertices(el)) cvc = convert.(Cartesian, vc) points = [Meshes.Point(Float64.(ustrip.((c.x, c.y, c.z)))...) for c in cvc] sel =Hexahedron(tuple(points...)) bb = Meshes.boundingbox(sel) # Collect voxels (cuboids) inside the hexahedron dbb = RegularGrid(bb.min, bb.max, dims = (10, 10, 10)) ***@***.*** nvertices(dbb) cuboids = [] for cl in nelements(dbb) c = element(dbb,cl) c ⊆ sel if c ⊆ sel # keep cuboid if inside element push!(cuboids, cl) end end — Reply to this email directly, view it on GitHub <#1260 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAZQW3PPVY6SSDD2UWVN6GT3Y65YXAVCNFSM6AAAAACJ4MCGZKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZTIMZTGQ2TQOBXGM> . You are receiving this because you were mentioned.Message ID: ***@***.***>