Instancing a custom number type and using in arrays

I’m trying to write FEM code for school and I’d like to use custom number types as a means of error-proofing. Essentially, I’d like to have a type difference between element ids, global node ids, local node ids, global dof ids, local dof ids, etc. Then arrays of these entries would be something like

elem_connectivity = Array{elem_id,2}(undef, num_local_nodes, num_elems)

etc. and it would be impossible to put the wrong kind of number in them. Additionally, I would like to make indexing the arrays require the right type of number, e.g

elem_connectivity[a::loc_node_id, e::elem_id] 

I’ve tried reading the manual on primitive types, but I seem to be unsure of how to begin creating instances of the types I create. For example

primitive type nodeID <: Signed 64 end
a = nodeID(100)

ERROR: MethodError: no method matching nodeID(::Int64)
Closest candidates are:
  nodeID(::T) where T<:Number at boot.jl:718
  nodeID(::Float16) where T<:Integer at float.jl:71
  nodeID(::Complex) where T<:Real at complex.jl:37
  ...
Stacktrace:
 [1] top-level scope at REPL[10]:1

Thoughts?

julia> primitive type nodeID 64 end

julia> nodeID(i::Integer) = reinterpret(nodeID, Int64(i))
nodeID

julia> nodeID(30)
nodeID(0x000000000000001e)

(If you’re not going to define arithmetic operations on it, don’t make it a subtype of Signed.)

But honestly, it sounds like your plan for a distinct type for each kind of index involves a lot of pain, little gain, and will result in non-idiomatic Julia code.

2 Likes

For your purposes, using struct types rather than primitive types makes the task easier and codes more concisely. To simplify your specifics (for clarity) there should be types ElementId, GlobalNodeId, LocalNodeId. Each of these should subtype Signed and otherwise behave as an Int.

This is elaborated code to help you see how to do this with structs; for example, you may not need abstract types beyond AbstractId.

abstract type AbstractId        <: Signed     end
abstract type AbstractNodeId    <: AbstractId end
abstract type AbstractElementId <: AbstractId end

struct ElementId    <: AbstractElementId
  id::Int
end

struct GlobalNodeId <: AbstractNodeId
  id::Int
end

struct LocalNodeId <: AbstractNodeId
  id::Int
end

# so values of these types display just as an Int does
for T in (:ElementId, :GlobalNodeId, LocalNodeId)
    @eval begin
         Base.show(io::IO, x::$T) = show(io, x.id)
    end
end

# to create concrete ids of each type, all valued `1`
element_id      = ElementId(1)
global_node_id  = GlobalNodeId(1)
local_node_id   = LocalNodeId(1)

You may have questions about what next / how to use these types to do _ .
Ask away. If you want to compare two ids of the same type those comparison operations need to be defined for each struct type.

import Base: <. <=, ==
for Op in (:(<). :(<=), :(==),)
    @eval begin
         $Op(x::ElementId, y::ElementId)       = $Op(x.id, y.id)
         $Op(x::GlobalNodeId, y::GlobalNodeId) = $Op(x.id, y.id)
         $Op(x::LocalNodeId, y::LocalNodeId)   = $Op(x.id, y.id)
    end
end
2 Likes

I’m ready for the pain - my advisor’s startup is developing an FEM code in C++ and what I described is how they’ve developed their code. I’m trying to follow this approach as much as possible in Julia. Sadly I don’t have access to their source so I don’t know exactly what they’re doing in C++. That being said, I’ll probably follow @JeffreySarnoff’s approach instead, because that seems a lot less painful.

@JeffreySarnoff I’m going to accept this as the solution to the question as posed, but I will probably take you up on your suggestion to ask away. I’m going to stew on / experiment with what you’ve described and will reply with any further questions.

It’s great that you have an existing FEM code to guide you in your choice of algorithms. However, realize that Julia and C++ are very different languages, and the choice of abstractions and code structure in C++ may not be a good template for a Julia implementation.

3 Likes