I hope I am asking the question in a correct way. I am working with JuAFEM and I have a example code given, which calculates the stiffness matrix elementwise.

@inbounds for q_point in 1:getnquadpoints(cellvalues_u) # loop over quadrature points
dΩ = getdetJdV(cellvalues_u, q_point)
for i in 1:n_basefuncs_u # loop over number of ansatz functions
BuT = symmetric(shape_gradient(cellvalues_u, q_point, i))
for j in 1:i
Bu = symmetric(shape_gradient(cellvalues_u, q_point, j))
Ke[i, j] += BuT ⊡ C ⊡ Bu * dΩ
end
end
end

Is there a way to calculate the shape gradient such that I do not need to calculate the stiffness matrix elementwise but for each element in matrix form. In particular I do not want to access each element of the stiffness matrix (Ke[i,j]) but calulate the matrix itself i.e.

@inbounds for q_point in 1:getnquadpoints(cellvalues_u) # loop GPi
dΩ = getdetJdV(cellvalues_u, q_point)
Bu = ? # Calculate somehow the shape gradient Bu and BuT by a function
BuT = transpose(Bu)
Ke += BuT * C * Bu * dΩ
end

# This is for 2D where symmetric second order tensors have 3 entries
Bu = zeros(3, n_basefuncs_u)
for i in 1:n_basefuncs_u #
Bu[:, i] = tovoigt(symmetric(shape_gradient(cellvalues_u, q_point, i)))
end

This notebook has an (old) example of three different ways of computing Ke in JuAFEM. I’m not sure if it still works.:

I don’t really understand why you would want to code like this though. What’s the advantage over the first way?

I have to agree with Kristoffer. In Julia there is no advantage to using big matrices to formulate the elementwise finite element expressions. I did that in https://github.com/PetrKryslUCSD/FinEtools.jl and I had to use buffers all over the place in order to avoid temporary allocations (which will eventually kill performance, if one is not careful).

JuAFEM does the right thing by calculating the element wise matrices entry by entry (most of the time, if I recall). StaticArrays really shine when used this way.