Hi,

I am confused by the segment function results:

using Meshes

s1 = Segment((1.5,1.5), (3.0,1.5))
s2 = Segment((1.5,1.5), (3.5,1.5))
s3 = Segment((3.0,1.0), (2.0,2.0))

s1 ∩ s3  # Point(2.4, 1.5), wrong
s2 ∩ s3  # Point(2.5, 1.5), correct

Is it a bug? Sorry if I made a mistake here; too tired to think...

Currently we are printing a big list with all vertices. We need to truncate the list like we do with other types.

I found an example where the intersection algorithm is failing:

julia> s1 = Segment((2.,1.), (1.,2.))
Segment{2,Float64}
  └─Point(2.0, 1.0)
  └─Point(1.0, 2.0)

julia> s2 = Segment((1.,0.), (1.,1.))
Segment{2,Float64}
  └─Point(1.0, 0.0)
  └─Point(1.0, 1.0)

julia> s1 ∩ s2
Point(1.0, 2.0)

I've read the paper multiple times, and yet I can't figure out what is wrong with the implementation. The algorithm implemented in intersecttype enters in case (II) as expected, but then it branches on the wrong intersection type at the last minute based on the non-zero winding number. I did the calculation by hand and all winding numbers are making sense. Also followed the steps in the VS Code debugger.

In page 332 of the paper, the author mentions the following condition inside case (II):

Can you spot any flaw in the implementation? Perhaps this is a typo in the paper?

Paper: https://sci-hub.se/10.1016/0167-8396(91)90019-8

cc: @serenity4

At some point, we will most likely want to compare meshing algorithms. For the moment, we may want to only evaluate the quality of triangulated meshes, to compare triangulation algorithms, but in the future other types of algorithms (e.g. quad meshing) may probably be implemented. Besides algorithmic time performance, there are many geometrical considerations that are of importance. This can be broken down into several characteristics. Here is a small general list to start with, which should be updated based on discussion and feedback:

• Mesh statistics:
  • Number of triangles (after triangulation, if considering non-triangle meshes).
  • Memory usage (for the total mesh; irrelevant between meshes with the same element types. Between triangle meshes with the same embedding dimension and coordinate type for example, it would be directly proportional to the number of elements.)
• Triangulation quality:
  • Sliverness (min, median, max)
• Application specific:
  • Number of poles
  • Number of involved edges in largest pole

In the scope of this package, this would be most relevant for testing, although I think some algorithms may directly benefit from this. External applications would definitely find this of interest.

Some geometric objects possess a natural parametric space. For example, segments linearly interpolate between two points, rays can be evaluated from the origin to the infinity, bezier curves can be evaluated between its endpoints... I had the idea of making these objects callable on their parametric inputs, so a segment could be called with t going from 0 to 1 to make a linear interpolation of its two points. Considering objects with a single parametric dimension it should be fine, unambiguous and easy to do. I can draft a PR to see how it goes for one-dimensional parametric objects.

When there is more than one parametric dimension, I am not sure how this should be implemented. For triangles, for example, it is possible to use barycentric coordinates (from 0 to 1), but also two-dimensional coordinates within a certain range (since this is a polygon after all, with parametric dimension 2). Feel free to comment if you have any ideas.

• Segment
• Ray
• BezierCurve
• Line
• Triangle
• Quadrangle

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

We have implemented the winding method of Balbes & Siegel 1990 for polygonal chains. It turns out that the method doesn't work with polygons that are not simple. I have added an empty implementation of Held, M. 1998 in chains.jl in case someone wants to give it a try. It should be a very simple calculation of triangular areas: https://link.springer.com/article/10.1007/s00453-001-0028-4

I would like to point out that MeshCore implements derivation of various adjacency relationships. I wonder if there is some synergy we could exploit here?

The Meshes.ears function returns the indices in a chain that correspond to ears as defined in the FIST paper. The current implementation assumes a CCW orientation of the chain and that is causing issues in the FIST triangulation.

We need to generalize the Meshes.ears implementation so that it works with both orientations.

I've migrated the trait coordinates! from GeoStatsBase.jl to Mehes.jl for computing the coordinates of the centroid of the i-th element of the mesh. Most algorithms follow this pattern:

function algorithm(mesh::{Dim,T})
  x = MVector{Dim,T} # pre-allocate memory
  for i in 1:nelements(mesh)
    coordinates!(x, mesh, i)
    # do something with x
  end
end

Now I realized that this pre-allocation of the coordinates with MVector doesn't make any difference in performance. I could have simply used a non-mutating version x = coordinates(mesh, i) in the loop for clarity. The benchmark bellow shows that there is no difference in CPU time, nor in memory allocations. Because the names coordinates! and coordinates are already taken, I am using the names centroid! and centroid in the benchmark:

using Meshes
using StaticArrays
using BenchmarkTools

function centroid!(x, Ω, i)
  for j in 1:embeddim(Ω)
    x[j] = 2.0*i + j
  end
  x
end

function centroid(Ω, i)
  SVector(ntuple(j->2.0*i+j, embeddim(Ω)))
end


function f1(Ω)
  Dim = embeddim(Ω)
  T = coordtype(Ω)
  x = MVector{Dim,T}(undef)
  s = zero(T)
  for i in 1:nelements(Ω)
    centroid!(x, Ω, i)
    s += sum(x) / Dim
  end
  s
end

function f2(Ω)
  Dim = embeddim(Ω)
  T = coordtype(Ω)
  s = zero(T)
  for i in 1:nelements(Ω)
    x = centroid(Ω, i)
    s += sum(x) / Dim
  end
  s
end

Ω = CartesianGrid(100,100)

@btime f1(Ω)
@btime f2(Ω)

8.738 μs (1 allocation: 16 bytes)
1.00025e8
8.932 μs (1 allocation: 16 bytes)
1.00025e8

In this proposal I would like to eliminate the mutating version coordinates! and replace the name coordinates by the name centroid(mesh, i) that allocates an SVector in the stack. This name is more appropriate and we can make this function return an actual Point instead of a raw vector.

I will prepare a PR with these modifications, but please let me know if you spot any failure in this analysis. After this PR is reviewed and merged, I will be able to release another series of versions for the GeoStats.jl stack to use the more explicit centroid name.

cc: @serenity4

As a quick workaround to include more features to the package, I've included a couple of rotation conventions from various software in geosciences under the conventions folder as we did with the GMSH ordering convention. These conventions should be replaced by actual types following the ReferenceFrameRotations.jl API so that we can easily convert back and forth between different rotation representations including quaternions, direction cosine matrices, etc.

I am super busy migrating functionality from GeoStatsBase.jl to here, so feel free to tackle this issue if you are free 👍🏽

The FIST algorithm is ready except for the recovery process that is needed when there exists self-overlapping edges.

Currently, we have UnstructuredMesh as the only data structure for unstructured meshes. It actually implements a Face-based structure that is useful for visualization. It is not useful at all for neighbor search or other kinds of queries.

We need to rename UnstructuredMesh to FaceMesh (Face-based mesh) or a similar name, and introduce other data structures such as Edge, HalfEdge, etc.

Currently we have plot recipes inside this repository as a temporary solution. In order to plot categorical values over meshes I had to add the CategoricalArrays.jl dependency temporarily. The plan is to move all the plot recipes to a separate repository with Makie.jl recipes instead. After that we can remove the CategoricalArrays.jl dependency from here.

I stumbled upon SciML's ColPrac which is a set of guidelines aimed towards package developers to facilitate collaborative work. I found it quite interesting and it includes a lot of good practices that we already seek to apply whenever possible (such as "PRs do one thing", "code changes should be tested", etc.). Maybe we should link to those (by adding the badge) in the README? That way, contributors will know more about the kind of expectations we have.

We can also choose to create our own document, usually done in the form of a CONTRIBUTING.md, which can either be self-contained or build upon ColPrac (or any other guide). Personally I'm fine using ColPrac alone, having in mind that these are guidelines and not hard requirements.

Hi,
I was experimenting with the library and saw that

julia> Quadrangle <: Polygon
false
julia> Triangle <: Polygon
false

which I found a bit odd at first. Then I realized that Polygon is a concrete type. I might be wrong, but wouldn't any 2-face be a polygon? I also understand that the current Polygon struct serves to construct any kind of polygon, simple or not, including triangles and quadrangles (geometrically, not Triangle nor Quadrangle instances).

Maybe would it be better to define a Polygon abstract type, and have the current Polygon concrete type be something like a PolygonArea or GenericPolygon?

I guess it could also be possible to use a new Polygon abstract type in a trait-based approach instead, which would apply to Polytope objects.

The Chain type is missing a few constructors. For example, we can't initialize one with a static vector:

julia> using Meshes, StaticArrays

julia> Chain(@SVector [Point2(0.,0.),Point2(1.,1.)])
ERROR: MethodError: no method matching Chain(::SVector{2, Point2})
Closest candidates are:
  (::Type{PL})(::TP...) where {PL<:Polytope, TP<:Tuple} at /home/belmant/.julia/dev/Meshes/src/polytopes.jl:29
  (::Type{PL})(::P...) where {PL<:Polytope, P<:Point} at /home/belmant/.julia/dev/Meshes/src/polytopes.jl:27
  Chain(::CircularArrays.CircularVector{T, A} where {T, A}) at /home/belmant/.julia/dev/Meshes/src/polytopes/chain.jl:15
  ...
Stacktrace:
 [1] top-level scope
   @ REPL[7]:1

nor can we with just points (as varargs), although the docstring mentions Chain(p1, p2, ..., pn) at the struct definition:

julia> Chain(Point2(0.,0.),Point2(1.,1.))
ERROR: MethodError: no method matching Chain(::SVector{2, Point2})
Closest candidates are:
  (::Type{PL})(::TP...) where {PL<:Polytope, TP<:Tuple} at /home/belmant/.julia/dev/Meshes/src/polytopes.jl:29
  (::Type{PL})(::P...) where {PL<:Polytope, P<:Point} at /home/belmant/.julia/dev/Meshes/src/polytopes.jl:27
  Chain(::CircularArrays.CircularVector{T, A} where {T, A}) at /home/belmant/.julia/dev/Meshes/src/polytopes/chain.jl:15
  ...
Stacktrace:
 [1] Chain(::Point2, ::Vararg{Point2, N} where N)
   @ Meshes ~/.julia/dev/Meshes/src/polytopes.jl:27
 [2] top-level scope
   @ REPL[8]:1

The second case should be fixed if the first one is too.

This is a reference to small inconsistencies that were left over from the #35 fix. Point constructors for all dimensions work now, but a small inconsistency remains with mixed-type inputs: Implicit conversion is disallowed for some constructors but not all:

julia> Point(0, 1.) # this is good
ERROR: MethodError: no method matching Point(::Int64, ::Float64)
Closest candidates are:
  Point(::V...) where V at REPL[7]:1
...

julia> Point((0, 1.)) # this is especially undesirable
Point{1, Tuple{Int64, Float64}}([(0, 1.0)])

julia> Point([0, 1.]) # implicit conversion
Point{2, Float64}([0.0, 1.0])

julia> Point{2,Float64}(0, 1.) # good
ERROR: MethodError: no method matching Point{2, Float64}(::Int64, ::Float64)
Closest candidates are:
  Point{Dim, T}(::V...) where {Dim, T, V} at REPL[5]:1
...

julia> Point{2,Float64}((0, 1.)) # needs a better error message
ERROR: DimensionMismatch("Can't construct a Point{2,Float64} with an input of length 1")
...

julia> Point{2,Float64}([0, 1.]) # implicit conversion
Point{2, Float64}([0.0, 1.0])

One should decide to either allow or forbid implicit conversion in all Point constructors and enforce that scheme. (But this is really not a pressing issue)

BTW: A few tests were created but disabled until consistency is enforced:

So far we never define a full polytope structure, but maybe it would be nice to have one. So that if you're interested in an object with all its vertices, edges, faces and so on (essentially, a mesh in 3D) without relying on conventions, there's a type available and you don't have to pass around all k-faces individually. We don't even need to encode the value of the vertices if we stick to an abstract polytope (which only describes connectivity, but not data).

I am considering adding the following trait for domain types in the package:

"""
    IsArrayLike{D}

Tells whether or not a domain type `D` behaves as an array of elements
in a geographic space. For example, a Cartesian grid can be seen as an
array of geometries (e.g. quadrangle, hexahedron). This information can
be used to optimize algorithms (e.g. slices, bounding boxes).
"""
@traitdef IsArrayLike{D<:Domain} <- isarraylike(D)
isarraylike(D) = false

Do you think this is an appropriate term for the concept? Basically I want to encode the fact that some domains like CartesianGrid are just a Cartesian arrangement of connected elements. So computing bounding boxes for example reduces to looking up the first and last elements instead of looping over all elements of the domain.

As discussed in #37, we could benefit from benchmarks for:

• informing the users of the costly and cheap functions
• tracking performance regressions during development

It can be implemented first as a simple benchmark folder, and possibly added as a GitHub action is desired.

As we write more algorithms, we can certainly benefit from more tests. We are happy to submit PRs.

cc: @henry2004y @PaulDebus

s1 = Segment(Point(4.0, 3.0), Point(0.0, -3.0))
s2 = Segment(Point(1.0, 2.0), Point(3.0, -2.0))

@code_warntype intersecttype(s1, s2)

returns (I leave only the beginning of the message)

julia> @code_warntype intersecttype(s1, s2)
MethodInstance for Meshes.intersecttype(::Segment{2, Float64, StaticArrays.SVector{2, Po
int2}}, ::Segment{2, Float64, StaticArrays.SVector{2, Point2}})
  from intersecttype(s1::Segment{2, T, V} where V<:AbstractArray{Point{2, T}, 1}, s2::Se
gment{2, T, V} where V<:AbstractArray{Point{2, T}, 1}) where T in Meshes at /home/skoffe
r/.julia/packages/Meshes/PDxIG/src/polytopes/segment.jl:69
Static Parameters
  T = Float64
Arguments
  #self#::Core.Const(Meshes.intersecttype)
  s1::Segment{2, Float64, StaticArrays.SVector{2, Point2}}
  s2::Segment{2, Float64, StaticArrays.SVector{2, Point2}}
Locals
  @_4::Union{Nothing, Tuple{Int64, Int64}}
  @_5::Union{Nothing, Tuple{Int64, Int64}}
  @_6::Tuple{StaticArrays.SOneTo{2}, Int64}
  @_7::Tuple{StaticArrays.SOneTo{2}, Int64}
  wyzero::Bool
  wxzero::Bool
  determinatey::Bool
  windingy::Float64
  determinatex::Bool
  windingx::Float64
  Q::Vector{Point2}
  y2::Point2
  y1::Point2
  x2::Point2
  x1::Point2
  ȳ::Point2
  x̄::Point2
  j@_21::Int64
  α::Float64
  β::Float64
  j@_24::Int64
  @_25::Bool
Body::Any

It looks like it originates from intersectcolinear, for two reasons:

1. this function has @error in one of the branches (and an output of @error has different type then all other branches).
2. It returns too many different types and Julia's compiler gives up on them.

Hi,

I found another potential bug for allocating a CartesianGrid object. Following the example in the document

a = CartesianGrid((-1.,),(1.,), dims=(100,))

This operation can never finish, no matter what values I input for dims. I also tried

a = CartesianGrid((100,),(10.,),(1.,))

Same issue I guess.

I stumbled upon a bug that occurs when a PolyArea has no inner chains (julia 1.6.0-rc2):

julia> vertices(poly)
ERROR: ArgumentError: reducing over an empty collection is not allowed
Stacktrace:
  [1] _empty_reduce_error()
    @ Base ./reduce.jl:299
  [2] mapreduce_empty(f::Function, op::Base.BottomRF{typeof(vcat)}, T::Type)
    @ Base ./reduce.jl:342
  [3] reduce_empty(op::Base.MappingRF{Meshes.var"#175#176", Base.BottomRF{typeof(vcat)}}, #unused#::Type{Tuple{Chain{2, Float32}, Any}})
    @ Base ./reduce.jl:329
  [4] reduce_empty_iter
    @ ./reduce.jl:355 [inlined]
  [5] reduce_empty_iter
    @ ./reduce.jl:354 [inlined]
  [6] foldl_impl
    @ ./reduce.jl:49 [inlined]
  [7] mapfoldl_impl
    @ ./reduce.jl:44 [inlined]
  [8] #mapfoldl#214
    @ ./reduce.jl:160 [inlined]
  [9] mapfoldl
    @ ./reduce.jl:160 [inlined]
 [10] #mapreduce#218
    @ ./reduce.jl:287 [inlined]
 [11] mapreduce
    @ ./reduce.jl:287 [inlined]
 [12] #reduce#220
    @ ./reduce.jl:456 [inlined]
 [13] reduce
    @ ./reduce.jl:456 [inlined]
 [14] vertices(p::PolyArea{2, Float32, Chain{2, Float32}})
    @ Meshes ~/.julia/dev/Meshes/src/polytopes/polyarea.jl:51
 [15] top-level scope
    @ REPL[12]:1

We should either add an empty vector as init argument to reduce, or use a different concatenation method for inner chains.

a1d408c is in violation of the MIT license under which the code from https://github.com/JuliaGeometry/GeometryBasics.jl is copied: see accompanying discussion in JuliaGeometry/GeometryBasics.jl#101 (comment).

Needing to proceed quickly with your own research is not an accepted reason for violating community standards and copyright notices -- please fix this issue immediately.

• To be clear: b3907dc (c.f. JuliaGeometry/GeometryBasics.jl#101 (comment)) is not the right way to fix it.
• Simon and Martjin have agreed that an acceptable way of resolving this issue is to fix the commit history across this repository to reflect that it is a fork of https://github.com/JuliaGeometry/GeometryBasics.jl.
• Failure to do so will result in an escalation of this issue to the Julia Stewards as a violation of https://julialang.org/community/standards/#give_credit.

The FIST implementation we have is reasonably robust, but it is still failing with some difficult polygonal areas. In our test suite we have a bunch of polygonal area examples in the data directory that are loaded for testing. Tests are passing for most polygonal areas in the family poly1...5.line, for example:

poly1.line

poly2.line

poly3.line

But some cases like poly3.line are not passing the visual test:

Notice how some edges cross. That shouldn't be the case.

With the latest changes, we can also handle polygonal areas with holes in the holes1...5.line family, for example:

hole2.line

hole1.line

But some cases like hole1.line have crossing edges too:

Notice how the left inner hole is crossed by edges.

Finally we have polygonal areas in the smooth1...5.line family that are working fine:

smooth2.line

smooth1.line

But similarly, some cases fail with crossing edges:

It would be nice to dive in debug mode to figure out what is happening. Both PolyArea and UnstructuredMesh can be plotted directly with our current plot recipes so that helps too in the debugging process.

If someone can read the FIST paper by Held and then read our implementation, which basically follows the exact same notation, that would be great. I think I've implemented all the ingredients except for one or two corner cases they discuss, which may be related to the issues we are seeing.

If someone wants to dive in, I can help explaining the current implementation in a video chat for example.

To my surprise the following does not error:

using Meshes: Segment

Segment((0, 0), (0, 1), (0, 2))

Is this intended?

There are a few geometric primitives already defined (ball, box, cylinder and sphere). Do we plan to complete this collection with other primitives (Bézier curves, splines, torus, circles...)? For example, I need to use Bézier curves, which are naturally defined from control points, should I make a PR to implement them here?

This is an open question. I would personally think that we may want to put everything in Meshes for now, and split parts of it later on when it becomes too big.

In various algorithms there is a pattern coordinates(centroid(domain, i)) that works well when the domain is a Mesh. However, when the domain is a PointSet the centroid isn't defined. We could simply add a method centroid(p::Point) = p to avoid over-engineering downstream. I am still thinking about the implications of doing this, but it sounds reasonable mathematically speaking?

Let me know if you feel that this doesn't make sense. I plan to add the method over the week as I find the time.

Hi @SimonDanisch , this is a follow up from yesterday's Discourse post where we shared that it would be nice to have an interface for general meshes, from simple Julia arrays (treated as regular meshes with zero storage for the coordinates) to general mixed-element meshes for various different purposes.

As an initial target, we could try to prove the concept that such interface would enable (1) geostatistical estimation/simulation of properties over meshes with GeoStats.jl, (2) physical simulation with available FEM solvers, and (3) plotting with Makie.jl

I understand that (3) is already done, and that perhaps (2) as well? Regarding (1), I would like to have your feedback on the interface copied/pasted below:

module Meshes

using StaticArrays: MVector

#--------------
# MESH TRAITS
#--------------
ndims(::Type{M}) where M = @error "not implemented"

ctype(::Type{M}) where M = @error "not implemented"

cbuff(m::Type{M}) where M = MVector{ndims(m),ctype(m)}(undef)

isstructured(::Type{M}) where M = Val{false}()

isregular(::Type{M}) where M = Val{false}()

COMPILE_TIME_TRAITS = [:ndims, :ctype, :cbuff, :isstructured, :isregular]

# default versions for mesh instances
for TRAIT in COMPILE_TIME_TRAITS
  @eval $TRAIT(::M) where M = $TRAIT(M)
end

elements(mesh::M) where M = @error "not implemented"

nelms(mesh::M) where M = length(elements(mesh))

#----------------------
# MESH ELEMENT TRAITS
#----------------------
coords!(x::AbstractVector, mesh::M, elm::E) where {M,E} =
  @error "not implemented"

function coords!(X::AbstractMatrix, mesh::M,
                 elms::AbstractVector) where M
  for j in 1:length(elms)
    coords!(view(X,:,j), mesh, elms[j])
  end
end

function coords(mesh::M, elm::E) where {M,E}
  x = cbuff(M)
  coords!(x, mesh, elm)
  x
end

function coords(mesh::M, elms::AbstractVector) where M
  X = Matrix{ctype(M)}(undef, ndims(M), length(elms))
  coords!(X, mesh, elms)
  X
end

coords(mesh::M) where M = coords(mesh, elements(mesh))

vertices(mesh::M, elm::E) where {M,E} = @error "not implemented"

nverts(mesh::M, elm::E) where {M,E} = length(vertices(mesh, elm))

#------------------------
# ELEMENT VERTEX TRAITS
#------------------------
coords!(x::AbstractVector, mesh::M, elm::E, vert::V) where {M,E,V} =
  @error "not implemented"

function coords!(X::AbstractMatrix, mesh::M, elm::E,
                 verts::AbstractVector) where {M,E}
  for j in 1:length(verts)
    coords!(view(X,:,j), mesh, elm, verts[j])
  end
end

function coords(mesh::M, elm::E, vert::V) where {M,E,V}
  x = cbuff(M)
  coords!(x, mesh, elm, vert)
  x
end

function coords(mesh::M, elm::E, verts::AbstractVector) where {M,E}
  X = Matrix{ctype(M)}(undef, ndims(M), nelms(mesh))
  coords!(X, mesh, elm, verts)
  X
end

end # module

It is a quite simple interface that covers all algorithms currently implemented in GeoStats.jl as far as I can tell. It consists of basic information about the dimension of the mesh and type of coordinates, an access function to an iterator of elements (e.g. elements could simply be the range 1:n), and functions to access the coordinates of the elements themselves (some pre-defined baricenter or centroid for different element types), and the coordinates of the vertices of the elements.

Does the current implementation in GeometryBasics.jl (and associated interface) cover this use case already? If not, can we work on it together to create a separate package (I suggest Meshes.jl) with a set of traits without implementation that covers the use case?

After that we could start extending this interface to include volume, surfacearea and other useful properties of elements that are constantly used in algorithms.

Keeping track of the fact that additional tests should be added for FIST, notably to test the correctness of the result. Additional tests introduced in #66 but eventually not merged in #68 included tests for the output regarding the number of triangles, connectivity coverage and vertex data.

Apologies if this is already possible and I just couldn't find it. It would be excellent to have methods for generating meshes from randomly ordered point clouds. Likewise, methods for refining them so they are "nice" (angles not too acute, etc.).

I've migrated some rotation functionality from GeoStatsBase.jl to Meshes.jl, and some of the tests are broken.

@rmcaixeta since you are much more familiar with this part of the codebase than I am, could you please take a look?

Here are the 3 broken tests:

I did my best to not modify the code, but had to split the rotation from the scaling functionality to gain more flexibility for the kinds of things we want to do in Meshes.jl.

Hi,

Glad to see this brand new package in the Julia community! Hopefully this will become the go-to choice for constructing a 2D/3D mesh in Julia!

I do have a question regarding parallelization: does it now support building a mesh in a distributed manner, or if not, will that be considered in the future?

As we discussed in #6, it would be nice to clearly separate what is related to FEM from polytope theory. Currently, a Face is used to denote the fact that a type is a polytope used in FEM. Then, common polytopes are defined as Face subtypes, with a convenience facets function that connects all polytope vertices into facets from the GMSH project conventions.

The specificity of Face subtypes (as it is in the code) that is related to FEM is just that a mapping is defined from the object type (triangle, tetrahedron...) to an array of Connectivity elements (holding each a set of points and the target polytope), according to the GMSH conventions.

A way to go would be to extract this FEM-specific knowledge and put it somewhere else, detached from the actual polytope. Having types for simple polytopes is useful because they have a specific structure that easily allows to build all k-faces from the vertices according to a convention. This convention should be specified by the user, with a default value that can be GMSH's.

Later on, a separate FEM part can define other polytopes that are only relevant in the field. For example, it can be some kinds of "triangles" with more than 3 vertices, possibly including one at the center. They are not triangles in a mathematical sense, because they have more than 3 edges; however, vertices are layed out so that, with a specific connectivity, edges can be defined that can be seen as a triangle. When a center point is defined, it is simply not connected to anything else.
That is, of course, if you consider all FEM points to represent actual vertices. They can just be "points of interest" that have nothing to do with the geometrical shape. In this case, an adequate sampling method, which rely on an ordering convention, should be defined for the relevant polytopes.

I just thought I'd let you know that the documentation pages are down. Perhaps the CI needs to be re-triggered?

We have a lot of functionality that is hidden in the source code. We need to add examples to illustrate some of the important features we already have.

Do you know what is happening here? How to fix this warning during precompilation?

julia> using Meshes
[ Info: Precompiling Meshes [eacbb407-ea5a-433e-ab97-5258b1ca43fa]
WARNING: Method definition _linear(D, Any) where {D} in module Meshes at /home/juliohm/.julia/packages/SimpleTraits/1wYi7/src/SimpleTraits.jl:331 overwritten on the same line (check for duplicate calls to `include`).
  ** incremental compilation may be fatally broken for this module **

Opening this issue to discuss the definitions and representations of lines and planes. The literature is vast on these objects and so we need to be careful to define something practical still. For some people, the term "line" refers to curved objects in non-Euclidean geometry. For others, a line is a straight line with parametric dimension 1. We need to list all these definitions we know, start giving them names, and finally agree on useful representations to implement.

@juliohm suggested perhaps it would be a good idea to put some of our geometry code from OpticSim.jl in Meshes.jl.

We currently support Bezier and NURBS patches, and Zernike, QType, and Chebyshev polynomials. I've never seen the last three surface types used anywhere except in optics so they may not be great candidates for inclusion in Meshes.jl. Bezier and NURBS patches are used in many fields. We implemented just enough functionality to make these two surface types useable in OpticSim. Not sure if that is enough for your purposes.

For this code to be useful to us in OpticSim we need partials with respect to the parameterizing variables. Our ray intersection returns intervals instead of intersection points. The intervals are used in our CSG computations. Would it make sense to move the partial and ray intersection code to Meshes.jl? It doesn't seem to fit with the philosophy of the rest of the package.

Have to also say that I'm 100% time committed between my job and maintaining OpticSim.jl so I wouldn't be able to maintain this code here on Meshes.jl. Is there a contributor here who wants to take on the job?

I've pushed a trait to the master branch under the name coordinates that returns the coordinates of the centroid of the i-th element of a mesh or point set:

julia> pset = PointSet((1,1), (2,2), (3,3))
3 PointSet{2,Int64}
  └─Point(1, 1)
  └─Point(2, 2)
  └─Point(3, 3)

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> coordinates(pset, 1)
2-element MVector{2, Int64} with indices SOneTo(2):
 1
 1

julia> coordinates(grid, 1)
2-element MVector{2, Float64} with indices SOneTo(2):
 0.5
 0.5

Do you think that the name coordinates is a fine name for the trait? The trait returns a raw vector as opposed to a Point object. Currently, we can't create mutable points, but perhaps we should? If we can create mutable points using a MVector as the internal coords field, then we could rename the function to centroid(grid, i) for example and compute the centroid with existing memory and as a Point.

It seems strange for Chain and Segment to be <: Polytope while BezierCurve <: Primitive. Do we want to separate these 'parametric primitives' into their own category?

Some tests are still using readpoly(Float64, fname) with an explicit Float64 type. This is because we first need to fix some issues in single precision. Opening this issue to keep track of this effort.

@serenity4 I was thinking about the name PolySurface I suggested last year but I think the name PolyArea is clearer. I will prepare a quick PR with the changes, ok?

Please let me know if you feel that the name PolyArea isn't clear.

Found an issue with the last branch of the intersection algorithm:

s1 = Segment(Point(0.94495744f0, 0.53224397f0), Point(0.94798386f0, 0.5344541f0))
s2 = Segment(Point(0.94798386f0, 0.5344541f0), Point(0.9472896f0, 0.5340202f0))

s1 ∩ s2
Segment{2,Float32}
  └─Point(0.9472896f0, 0.5340202f0)
  └─Point(0.94798386f0, 0.5344541f0)

In this plot s1 is the black segment and s2 is the red segment. The intersection algorithm reaches the CASE (III) correctly, and then tries to decide whether or not these two segments are co-linear:

elseif !determinatex && !determinatey # CASE (III)
    if !isapprox((x2 - x1) × (y2 - y1), zero(T), atol=atol(T))
      # configuration (3)
      CornerTouchingSegments(intersectpoint(s1, s2))
    else
      # configuration (3), (4) or (5)
      intersectcolinear(s1, s2)
    end
  end

I am using the cross product to check for co-linearity, but the tolerance seems to be too tight. The atol(T) should return a reasonable tolerance for distances, not areas or volumes, I think that is causing the issue. Based on the plot, we would like to return configuration (3) because the segments aren't overlapping and aren't co-linear.

How should we set the tolerance in this case? Or should we use a different condition to check for co-linearity that doesn't use the cross product?

Hello,

When playing with the code, I got a question about the type system in the package. For example,

using Meshes

points = Point2[[1.0,1.0],[3.0,1.0]]
a = Segment(points[1], points[2])
typeof(a) == Segment # false?
typeof(points[1]) == Point # false?
typeof(points[1]) == Point2 # true?

Do I need to specify all the parameters for the type, e.g. typeof(a) == Segment{2,Float64,StaticArrays.SArray{Tuple{2},Point{2,Float64},1,2}}? Why is Point not a general supertype for all points?

After installing v0.3, using Meshes produces the following error:

ERROR: ArgumentError: Package Meshes [eacbb407-ea5a-433e-ab97-5258b1ca43fa] is required but does not seem to be installed:
 - Run `Pkg.instantiate()` to install all recorded dependencies.

Stacktrace:
 [1] _require(::Base.PkgId) at .\loading.jl:999
 [2] require(::Base.PkgId) at .\loading.jl:928
 [3] require(::Module, ::Symbol) at .\loading.jl:923

I also tried Pkg.dev("Meshes"); using Meshes, which failed a different way:

[ Info: Precompiling Meshes [eacbb407-ea5a-433e-ab97-5258b1ca43fa]
ERROR: LoadError: LoadError: cannot replace module Meshes during compilation
Stacktrace:
 [1] include(::Function, ::Module, ::String) at .\Base.jl:380
 [2] include at .\Base.jl:368 [inlined]
 [3] include(::String) at C:\Users\...\.julia\dev\Meshes\src\Meshes.jl:5
 [4] top-level scope at C:\Users\...\.julia\dev\Meshes\src\Meshes.jl:24
 [5] include(::Function, ::Module, ::String) at .\Base.jl:380
 [6] include(::Module, ::String) at .\Base.jl:368
 [7] top-level scope at none:2
 [8] eval at .\boot.jl:331 [inlined]
 [9] eval(::Expr) at .\client.jl:467
 [10] top-level scope at .\none:3
in expression starting at C:\Users\...\.julia\dev\Meshes\src\meshes.jl:5
in expression starting at C:\Users\...\.julia\dev\Meshes\src\Meshes.jl:24
ERROR: Failed to precompile Meshes [eacbb407-ea5a-433e-ab97-5258b1ca43fa] to C:\Users\...\.julia\compiled\v1.5\Meshes\FuRcu_a2hQM.ji.
Stacktrace:
 [1] error(::String) at .\error.jl:33
 [2] compilecache(::Base.PkgId, ::String) at .\loading.jl:1290
 [3] _require(::Base.PkgId) at .\loading.jl:1030
 [4] require(::Base.PkgId) at .\loading.jl:928
 [5] require(::Module, ::Symbol) at .\loading.jl:923

This also happens when calling the fully-specified constructor Point{1,Int}(SA[1]). (I don't think there is right now any way to define a 1-dimensional point!). The simplest way to fix this is probably to define an inner constructor.