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