wip on using Tensors

[KristofferC]

Still quite WIP because there is quite a lot of churn to go through when changing this.

I believe something like this is a good way forward. It will take a while to push through but in the end it should be worth it.

Changing

function jacobian(B, X, xi)
    dB = eval_dbasis(B, xi)
    dim1, nbasis = size(B)
    dim2 = length(first(X))
    J = zeros(dim1, dim2)
    for i=1:dim1
        for j=1:dim2
            for k=1:nbasis
                J[i,j] += dB[i,k]*X[k][j]
            end
        end
    end
end

to

function jacobian(B::AbstractBasis, X::Vector, xi)
    dB = eval_dbasis(B, xi)
    J = zero(Tensor{2, dim(B)})
    for i = 1:length(X)
        J += dB[i] ⊗ X[i]
    end
    return J
end

is a pretty good example I feel.

[TeroFrondelius]

Yes definitely. I talked with @ahojukka5 today and we believe there is a big potential performance gain in this approach.

[KristofferC]

For multithreading we really need to reduce allocations to a minimum. When Julia's GC run, all threads stop so allocating is extra bad from a multithreading point of view. So starting to use these type of static data structures will help with that I feel.

[TeroFrondelius]

We originally made a decision that we should stick in ascii characters. I think now is the correct time to review our decision. Here the special operator syntax really helps readability.

[KristofferC]

I can change it to otimes(dB[i], X[i]) but I personally kind of like it. As long as it isn't forced upon users it can perhaps be OK.

[ahojukka5]

I like otimes more.

[KristofferC]

Ok, will change to that then.

And instead of

ϵ_i ⊡ C ⊡ ϵ_j

we do

dcontract(epsilon_i, dconstract(C, epsilon_j))

I guess.

[ahojukka5]

The problem with utf-8 characters, in general, is that people do not know how to make them. I did find someone using nabla in her code, and when I asked how did she do that character, the answer was that she copypasted it everytime from somewhere else (did not know about \nabla + tab). My opinion is that utf-8 characters look cool but there are more downsides than benefits when using them.

[KristofferC]

I updated this to use Tensors throughout and I rewrote some things just how I would have done it. It won't be directly usable with the rest of JuliaFEM but perhaps it can be used for some inspiration in using Tensors.jl in the future.

As a performance comparison:

Before

b = BasisInfo(Quad4)
X = ((0.0,0.0), (1.0,0.0), (1.0,1.0), (0.0,1.0))
xi = (0.0, 0.0)
@btime eval_basis!($b, $X, $xi)
  203.700 ns (0 allocations: 0 bytes)

After

b = BasisInfo(Quad4)
X = Vec.([(0.0,0.0), (1.0,0.0), (1.0,1.0), (0.0,1.0)])
xi = Vec(0.0, 0.0)
@btime eval_basis!($b, $X, $xi)
  38.464 ns (0 allocations: 0 bytes)

I still haven't integrated the manifold + curve thing for the Jacobian since I am a bit unexperienced with that.

[ahojukka5]

Would it be more Julia-style to use macro instead of create_and_eval?

[KristofferC]

Maybe, you aren't really doing any syntax transformations here though so perhaps not? I dont think it matters so much, I like to avoid macros as much as possible though since you cant really know what they do unless you carefully read the docs for them.

[KristofferC]

Added a few more performance tweaks, shaved of 10 ns.

[ahojukka5]

For 2d surfaces in 3d space we need Jacobian to determine the normal direction. (I'm not sure is this terminology right.) Then detJ is the size factor needed for numerical integration.

In practice:

using FEMBase, FEMBasis
X = Dict(1=>[0.0,0.0,0.0], 2=>[1.0,0.0,1.0], 3=>[1.0,1.0,1.0], 4=>[0.0,1.0,0.0])
element = Element(Quad4, (1, 2, 3, 4))
update!(element, "geometry", X)
bi = BasisInfo(Quad4)
element_info!(bi, element, (0.0,0.0), 0.0)

As a result, the cross product of the column vectors of the Jacobian is a normal direction and "size factor" is the norm of the normal vector:

t1 = bi.J[1,:]
t2 = bi.J[2,:]
n = cross(t1, t2)/norm(cross(t1,t2))

# output

3-element Array{Float64,1}:
 -0.7071067811865475
  0.0
  0.7071067811865475

detJ1 = norm(cross(n1,n2))
detJ2 = sqrt(det(bi.J*bi.J'))
detJ1 == detJ2 == bi.detJ == sqrt(2)/4

#output

true

For curves,

So is this the r'(t) in here?

Some links:

• https://math.stackexchange.com/questions/1246990/integrating-over-an-embedded-manifold-jacobian-factor
• https://en.wikipedia.org/wiki/Surface_integral
• https://en.wikipedia.org/wiki/Line_integral

[ahojukka5]

Should we consider implementing also functions to return a single shape function, given its number? get_basis(xi, 1) and so on?

[coveralls]

Coverage decreased (-41.08%) to 58.257% when pulling 3601116 on kc/tensors into 02da5c5 on master.