I agree with @Mauro on the definition of `FEMesh`.
When @Kristo writes `elements::Vector{FElement}`
you implied that `elements` can be an heterogeneous vector, which contains
`LinTrid` and `LinQuad` at the same time.
So the compiler has to test each component of the vector, which makes the
optimization of Julia totally useless.
That is why we need *parametric type*.
On Saturday, March 14, 2015 at 11:06:06 PM UTC+1, Mauro wrote:
>
> > type LinTrig <: FElement
> > v1::Int
> > end
>
> Here you should probably use `immutable`, that way it will be stored as
> efficiently as ints.
>
> > type LinQuad <: FElement
> > v1::Int
> > v2::Int
> > end
> >
> >
> > type Node
> > c::Float64
> > end
> >
> >
> > type FEMesh
> > node::Vector{Node}
> > elements::Vector{FElement}
> > end
>
> This is type-unstable and will be slow. For a single type of element
> this type should be written as:
>
> type FEMesh{T<:FElement}
> node::Vector{Node}
> elements::Vector{T}
> end
>
> As you suggested, I also think, to have several element types, splitting
> up the mesh into separate parts would be best.
>
> > function assembleK(mesh::FEMesh)
> > for element in mesh.elements
> > println(element.v1)
> > end
> > end
> >
> > Now, running @code_warntype gives
> >
> > @code_warntype(assembleK(mesh))
> >
> > Variables:
> > mesh::FEMesh
> > #s2::Int64
> > element::FElement
> >
> > Body:
> > begin # none, line 2:
> > GenSym(0) =
> (top(getfield))(mesh::FEMesh,:elements)::Array{FElement,1}
> > #s2 = 1
> > GenSym(2) = #s2::Int64
> > GenSym(4) = (top(arraylen))(GenSym(0))::Int64
> > unless
> >
> (top(box))(Bool,(top(not_int))((top(slt_int))(GenSym(4),GenSym(2))::Bool))
> > goto 1
> > 2:
> > GenSym(10) = (top(arrayref))(GenSym(0),#s2::Int64)::FElement
> > GenSym(11) = (top(box))(Int64,(top(add_int))(#s2::Int64,1))
> > element = GenSym(10)
> > #s2 = GenSym(11) # line 3:
> > println((top(getfield))(element::FElement,:v1)::Any)::Any
> > 3:
> > GenSym(6) = #s2::Int64
> > GenSym(8) = (top(arraylen))(GenSym(0))::Int64
> > unless
> >
> (top(box))(Bool,(top(not_int))((top(box))(Bool,(top(not_int))((top(slt_int))(GenSym(8),GenSym(6))::Bool))))
>
>
> > goto 2
> > 1:
> > 0:
> > return
> > end::Void
> >
> >
> > So I have a type instability that will cause a lot of allocations and
> slow
> > the code down. If I instead stick to one element type and instead define
> > mesh to contain a vector of the concrete type like this
> >
> > type FEMesh
> > node::Vector{Node}
> > elements::Vector{LinTrig}
> > end
> >
> > then I don't have an instability any more and my code runs much faster.
> >
> > My questions is basically how I would easiest remove the instability
> > without "hard coding" the type the elements will have.
> >
> >
> > On Saturday, March 14, 2015 at 10:35:04 PM UTC+1, Sisyphuss wrote:
> >>
> >> Hello Kristoffer, I don't quite understand what you are doing. But by
> >> judging from the title of your post, I'd like to share with you my
> opinion,
> >> and I hope it would be helpful.
> >>
> >> The (multiple) dispatch used by Julia is a dynamic one, that is, the
> >> dispatch is preformed at run time instead of compiling time (or "code
> write
> >> time" as you mentioned). So you shouldn't have any worry about that, as
> >> long as your callee function can handle different concrete types. If it
> is
> >> not the case or it can't be achieved in a simple way, it is advised to
> >> write separate functions (a.k.a. generic function) for different
> concrete
> >> types.
> >>
> >>
> >>
> >> On Saturday, March 14, 2015 at 9:28:48 PM UTC+1, Kristoffer Carlsson
> wrote:
> >>>
> >>> Hello,
> >>>
> >>> I am writing some code for Finite Elements in Julia. As usual, I have
> a
> >>> mesh with elements and nodes. However, I do not want to restrict
> myself to
> >>> only using one type of element per mesh and therefore I have a
> hierarchy of
> >>> element types along the lines of:
> >>>
> >>> abstract FElement
> >>>
> >>> type LinTrig <:FElement end
> >>>
> >>> type LinQuad <: FElement end
> >>>
> >>> etc
> >>>
> >>> My mesh then consists of a Vector of FElements.
> >>>
> >>> Since I have a Vector of an abstract type I get a performance hit when
> I,
> >>> for example, assemble the global stiffness matrix. The compiler does
> not
> >>> know what type of element it will find and have to check it in each
> >>> iteration. For simple material models the overhead from the type
> checking
> >>> is 10 times larger than calculating the elements stiffness matrix
> >>> contribution.
> >>>
> >>> Now, when I am writing the code for the assembler I can't know what
> >>> elements that the mesh will contain. However, *at the actual call* to
> >>> the assembly function I know what element types I have in the mesh and
> that
> >>> these will never change. I should then be able to give the compiler
> that
> >>> information and have it generate efficient code.
> >>>
> >>> I mean, I can do this manually. Just check for the element types,
> >>> generate one vector of elements for each element type and call the
> assembly
> >>> function for each separate vector of elements. All of these types will
> be
> >>> known.
> >>>
> >>> Is this possible to do somehow?
> >>>
> >>> I hope I made sense.
> >>>
> >>> Best regards,
> >>> Kristoffer Carlsson
> >>>
> >>
>
>
