Is this an ok place to ask questions? Not sure if there is a gitter channel

[goedman]

Hi Kristoffer,

I'm working on 'translating' a FEM toolkit (http://ptfem.github.io/PtFEM.jl/latest/INTRO.html) that I have been using in Fortran & R for quite some years (these days mainly for directional drilling applications for geothermal holes and occasionally buckling). The book contains 7 chapters demonstrating a variety of applications of the PtFEM toolkit.

I've converted about 2 1/3 chapters out of the 7 to Julia and started (in parallel) to explore the possibility of making JuaFEM the core layer, i.e. underneath what the book describes as the program layer (p41.jl, p52.jl, etc) and complete that item maybe earlier then I had anticipated in my TODO list.

Not sure if the exploration shows it is possible and if so if going that route is a better approach, but looking at JuAFEM evaluating this option is certainly worth a 1 or 2 week effort from my side.

Apologies for the long intro. Finally, my questions:

1. Exploring JuAFEM's grid.jl, a fairly simple first question: Is in JuAFEM a Cell always synonymous with a finite element? In the latest docs under "Implementation in JuAFEM" it does state "elements or cells". PtFEM distinguishes structural elements like in your examples the cantilever beam.

2. Looking at interpolations.jl, the shape functions look similar (they are numbered differently). Is the intention that a package like PtFEM.jl add interpolations like the 20-node hexahedron or would you prefer PRs?

Still early on, but I sure like what I'm seeing in JuAFEM.jl!

Regards, Rob

[fredrikekre]

Hello and welcome!

[...] explore the possibility of making JuaFEM the core layer

I think that JuAFEM is perhaps more of a backend/core layer than a complete finite element software, and that is also the intention. So I definitely think it will work well as a core layer, and I am sure there are things that can be improved to make that simpler as well.

1. Exploring JuAFEM's grid.jl, a fairly simple first question: Is in JuAFEM a Cell always synonymous with a finite element? In the latest docs under "Implementation in JuAFEM" it does state "elements or cells". PtFEM distinguishes structural elements like in your examples the cantilever beam.

The idea is that a Cell is just the geometric description (line, triangle, etc) defined only by the nodal coordinates. This means that a Cell does not contain any other information (like degrees of freedom). From a quick look at your code, StructuralElements seem to be collections of finite elements, for instance a Rod, which essentially is built up by several simple line elements? If I got that right, then we don't have any structural elements in JuAFEM. It shouldn't be hard to build your structural elements based on this though.

2. Looking at interpolations.jl, the shape functions look similar (they are numbered differently). Is the intention that a package like PtFEM.jl add interpolations like the 20-node hexahedron or would you prefer PRs?

It would be great to support 20-node hexahedron! Feel free to open a PR.

[goedman]

Thanks a lot Fredrik, that helps a lot. Correct, StructuralElements are a mesh/grid of FiniteElements, as grid-generator() does in the cantilever example. PtFEM.jl's FiniteElements do hold the DOFs per node. I'll study JuAFEM a bit further and construct some test cases - just to make sure I get it - before opening a PR.

[mohamed82008]

Hi,

It seems like this issue can be used for asking design questions if no one here minds, otherwise I can move it to another issue or something. I wonder why in the following piece of code starting at line 90 in FEValues/cell_values.jl that you used a Vec/first order tensor to represent shape functions and a second order tensor to represent the Jacobian, when in fact you are using one of the entries only (the others are 0) in the shape function first order tensor, and one row of the Jacobian matrix from which the second order tensor is made? Was this just to pave the road to using vector shape functions and scalar dofs at each point, which is often used in non-mechanical FEA, as opposed to scalar shape functions and vector dofs which are the standard in mechanical FEA? Or is there any other reason I don't foresee? Thanks in advance! Great package btw and very educational to read.

# CellVectorValues
struct CellVectorValues{dim, T <: Real, refshape <: AbstractRefShape, M} <: CellValues{dim, T, refshape}
    N::Matrix{Vec{dim, T}}
    dNdx::Matrix{Tensor{2, dim, T, M}}
    dNdξ::Matrix{Tensor{2, dim, T, M}}
    detJdV::Vector{T}
    M::Matrix{T}
    dMdξ::Matrix{Vec{dim, T}}
    qr_weights::Vector{T}
end

function CellVectorValues(quad_rule::QuadratureRule, func_interpol::Interpolation, geom_interpol::Interpolation=func_interpol)
    CellVectorValues(Float64, quad_rule, func_interpol, geom_interpol)
end

function CellVectorValues(::Type{T}, quad_rule::QuadratureRule{dim, shape}, func_interpol::Interpolation,
        geom_interpol::Interpolation=func_interpol) where {dim, T, shape <: AbstractRefShape}

    @assert getdim(func_interpol) == getdim(geom_interpol)
    @assert getrefshape(func_interpol) == getrefshape(geom_interpol) == shape
    n_qpoints = length(getweights(quad_rule))

    # Function interpolation
    n_func_basefuncs = getnbasefunctions(func_interpol) * dim
    N    = fill(zero(Vec{dim, T})       * T(NaN), n_func_basefuncs, n_qpoints)
    dNdx = fill(zero(Tensor{2, dim, T}) * T(NaN), n_func_basefuncs, n_qpoints)
    dNdξ = fill(zero(Tensor{2, dim, T}) * T(NaN), n_func_basefuncs, n_qpoints)

    # Geometry interpolation
    n_geom_basefuncs = getnbasefunctions(geom_interpol)
    M    = fill(zero(T)            * T(NaN), n_geom_basefuncs, n_qpoints)
    dMdξ = fill(zero(Vec{dim, T})  * T(NaN), n_geom_basefuncs, n_qpoints)

    for (qp, ξ) in enumerate(quad_rule.points)
        basefunc_count = 1
        for basefunc in 1:getnbasefunctions(func_interpol)
            dNdξ_temp, N_temp = gradient(ξ -> value(func_interpol, basefunc, ξ), ξ, :all)
            for comp in 1:dim
                N_comp = zeros(T, dim)
                N_comp[comp] = N_temp
                N[basefunc_count, qp] = Vec{dim, T}((N_comp...))

                dN_comp = zeros(T, dim, dim)
                dN_comp[comp, :] = dNdξ_temp
                dNdξ[basefunc_count, qp] = Tensor{2, dim, T}((dN_comp...))
                basefunc_count += 1
            end
        end
        for basefunc in 1:n_geom_basefuncs
            dMdξ[basefunc, qp], M[basefunc, qp] = gradient(ξ -> value(geom_interpol, basefunc, ξ), ξ, :all)
        end
    end

    detJdV = fill(T(NaN), n_qpoints)
    MM = Tensors.n_components(Tensors.get_base(eltype(dNdx)))

    CellVectorValues{dim, T, shape, MM}(N, dNdx, dNdξ, detJdV, M, dMdξ, quad_rule.weights)
end

[fredrikekre]

Was this just to pave the road to using vector shape functions and scalar dofs at each point

Yes. Essentially, for a vector valued problem you have two options:

1. Use CellScalarValues (scalar shape functions) in combination with "vector dofs"
2. Use CellVectorValues (vector valued shape functions) in combination with "scalar dofs"

[mohamed82008]

But this is not usable yet I assume, because only one shape function is non-zero?

Edit: It actually seems fine.

[fredrikekre]

But this is not usable yet I assume

It is usable. Option 2. in my last post is used here: https://github.com/KristofferC/JuAFEM.jl/blob/master/examples/cantilever.ipynb

[mohamed82008]

Yup I realized it was ok after writing the comment. Thanks.

[mohamed82008]

Can you please provide a reference for this part of the code in line 49 in FEValues/face_integrals.jl and the similar parts there? This seems to be computing the normal vector to a face weighted by the area of the face using the whole element's geometric Jacobian matrix and the face index. Thank you.

function weighted_normal(J::Tensor{2, 2}, ::FaceValues{2, T, RefCube}, face::Int) where {T}
    @inbounds begin
        face == 1 && return Vec{2}(( J[2,1], -J[1,1]))
        face == 2 && return Vec{2}(( J[2,2], -J[1,2]))
        face == 3 && return Vec{2}((-J[2,1], J[1,1]))
        face == 4 && return Vec{2}((-J[2,2], J[1,2]))
    end
    error("unknown face number: $face")
end

[KristofferC]

See e.g. https://scicomp.stackexchange.com/a/26889

[mohamed82008]

Thanks, I will check it out tomorrow!

[mohamed82008]

So I just made it out of the FEValues folder, which is probably the Mathiest part of the whole package, and I saw you had some plans for vector shape functions and vector dofs in the common_values.jl file, making function evaluation, gradient and divergence functions for them. So what is on the radar next for this project? Is there any extensions or roadmap or anything like that? I still have some way to read the whole package and be capable of contributing but I am just curious to know what's on the table.

My main application is topology optimization so I will probably have a lot on my own table but I would love to contribute here as well if I can, and if there is something of common interest.

[KristofferC]

If you use the VectorValues type of the FEValues then you get out vectors from e.g. shape_value and a second order tensor from shape_gradient. This is used in the Hyperelasticity example https://github.com/KristofferC/JuAFEM.jl/blob/master/examples/hyperelasticity.ipynb.

One thing to work on is to add an example of topology optimization to the notebooks. Perhaps you are interested in https://github.com/KristofferC/MMA.jl by the way. I haven't looked at it for a long while though so it probably needs updating.

As to roadmap, the current main things that are planned are the following:

• Fix up https://github.com/KristofferC/JuAFEM.jl/pull/138. The implementation for this is quite complicated so it is nothing that I would recommend to work on, but it would allow us to move away from using isoparametric elements.

• Look over the Grid datastructures and see if they can be made better / more flexible.

• Perhaps add support for hanging nodes. This would make implementing adaptive mesh refinement easier.

[mohamed82008]

I found MMA.jl and it was next on my reading list, thanks for pointing it out. I think starting with the topology optimization examples would be great. I like the grid improvements as well because topology optimization could benefit from adaptive meshing, and my work on https://github.com/mohamed82008/VTKDataTypes.jl will be relevant here. Let's see how much I can do.

[mohamed82008]

In line 99 in Grid/grid.jl, the following function is defined as:

getfacelist(grid::Grid) = getfacelist(eltype(grid.cells))

but I cannot find the definition of getfacelist for any Type{<: Cell}. Am I missing something?

[KristofferC]

The Grid stuff has been a bit in flux. That method is probably unused. There are methods for Interpolations

[mohamed82008]

Ya I found the ones for interpolations but not this one.

[mohamed82008]

Is there any specific reason for using N and M in the type definition of Cell? What do you think of adding a reference_shape and order as parameters in the type Cell instead of N and M which can be inlined functions of the types, keeping dim? Then a Cell subtype can be used as a parameter in the Interpolation type together with the interpolation order. Then getfacelist which is really a function of the Cell type can be defined for each Cell type and for an Interpolation using its Cell parameter. Also the Serendipity interpolation may either be associated with a full cell used in Lagrange interpolation (and in VTK), with all its nodes present, or a new Cell type can be made with fewer nodes specifically for the Serendipity to save on memory. Currently, grids and interpolations seem like separate islands but they are really just using different encodings of similar things.

Thinking out loud (well maybe not too loud): for the non-isoparametric case, the order of the cell must be greater than or equal to the order of the geometric interpolation and function interpolation. This means that we could provide a way to view a higher order cell as a lower order cell in case the geometric interpolation was lower order and only needed to get corner points to evaluate the function or Jacobian, and another function to interpolate the dofs for the extra nodes in the higher order cell in case the function interpolation was of lower order than the geometric interpolation. The quadrature can remain unchanged as it is not really dependent on the cell type but only on the reference shape, because it is used to approximate integrals over the reference shape regardless of the function inside the integral.

Any comments on all the above? Is this making it more organized or more complicated?

[mohamed82008]

In lines 205 and 206 in Dofs/DofHandler.jl, shown below, why is the second line useful? If I understand correctly, you want N to be an upper bound on the number of non-zero values in the matrix to allocate memory, so is there a case where the first line is not big enough? Seems sufficient to me.

N = sym ? div(n*(n+1), 2) * ncells : n^2 * ncells
N += ndofs(dh) # always add the diagonal elements

[fredrikekre]

If I understand correctly, you want N to be an upper bound on the number of non-zero values in the matrix to allocate memory

Not quite, N is the total number of elements that we will push to I and J. So N = sym ? div(n*(n+1), 2) * ncells : n^2 * ncells we add the diagonal covers https://github.com/KristofferC/JuAFEM.jl/blob/7e4aea5b8e3217a119fb3ac0cb9ce7b29ce63546/src/Dofs/DofHandler.jl#L216-L217 and N += ndofs(dh) covers https://github.com/KristofferC/JuAFEM.jl/blob/7e4aea5b8e3217a119fb3ac0cb9ce7b29ce63546/src/Dofs/DofHandler.jl#L221-L222

[mohamed82008]

Hmm, I see. So for the code below, the first loop visits every entry in the stiffness matrix and records the row and column indices. Some of these will be repetitive but the function sparse will take care of that. But why is the second loop necessary? Shouldn't all the diagonal entries be visited in the first loop already?

    for element_id in 1:ncells
        celldofs!(global_dofs, dh, element_id)
        @inbounds for j in 1:n, i in 1:n
            dofi = global_dofs[i]
            dofj = global_dofs[j]
            sym && (dofi > dofj && continue)
            push!(I, dofi)
            push!(J, dofj)
        end
    end
    for d in 1:ndofs(dh)
        push!(I, d)
        push!(J, d)
    end
    V = zeros(length(I))
    K = sparse(I, J, V)

[KristofferC]

It is just to be 100% sure that we have stored entries on the diagonal since we use that fact later.

[mohamed82008]

Actually if the first loop skips even 1 diagonal element then the stiffness matrix can never be positive definite, which is a contradiction.

[KristofferC]

Contradiction of what?

[mohamed82008]

Logical contradiction, because it is known/proven that the stiffness matrix is positive definite.

[KristofferC]

I guess you haven't used mixed elements or Lagrange multipliers? If my problems where always positive definite, I would be so happy.

[mohamed82008]

Oh ok then it's my lack of knowledge then sorry.

[KristofferC]

As an example,

[mohamed82008]

Great I know about contact constraints and Lagrangian multipliers but I guess not enough to realize that K will be positive semidefinite. Btw the whole matrix above is indefinite, so they mean K is positive semidefinite I assume not the whole thing?

[mohamed82008]

Actually the only things that can make K not positive definite is if the stiffness tensor is a bit funny, having 0 or even negative values in some directions or so, or if some nodes are not part of any element. The latter would make it positive semidefinite. So how did I know the above matrix is indefinite just by looking? Referring to 4.2.8 in the 3rd edition of "Matrix Computations" by Gene H. Golub and Charles F. Van Loan, if the matrix is positive semidefinite and any diagonal element is 0, then all the elements of the row or column of that diagonal must also be 0. The rationale of the proof is as follows. If K is positive (semi-)definite, then it is nothing but the Gram matrix of the rows of the Cholesky factor. That is K_ij = <v_i,v_j> where v_i is the vector made from the ith row of the Cholesky factor. So if the ith diagonal element is 0, then <v_i,v_i> is 0, so v_i is 0, so <v_i, v_j> are 0 for all j.

The same argument can be applied on minus the matrix, so the above matrix is also not negative definite and not negative semidefinite and it can only be indefinite. The only case for it to possibly be positive semidefinite is if A were 0 which it isn't.

Now this is getting severely off topic so you can keep the loop if you want!

[natschil]

Hi, if I understood it correctly, it is ok to ask general questions about JuAFEM here, so I hope no-one will mind if I add this one:

What is the role of the geometric interpolation functions? I'm trying to understand the following code from FEValues/cell_values:

      fecv_J = zero(Tensor{2, dim})
        for j in 1:n_geom_basefuncs
            fecv_J += x[j] ⊗ cv.dMdξ[j, i]
        end

As far as I can tell, x[j] is the j-th vertex of the cell, and cv.dMdξ[j, i] is the gradient of the j-th geometric interpolation function at the current quadrature point. Could someone be so kind as to point me to a reference that describes what the geometric interpolation functions are, and why this tensor is used for the transformation of the volume?

Best Regards, Nathanael

[KristofferC]

The geometric interpolation is used to compute for example the Jacobian J such that det(J)is the volume element that is needed in the quadrature integration. 17.2.1 here https://www.colorado.edu/engineering/CAS/courses.d/IFEM.d/IFEM.Ch17.d/IFEM.Ch17.pdf shows a bit about it.

The geometric interpolation is also used if you want to compute global coordinates in e.g. a quadrature point.

[KristofferC]

For further questions, can we use the #fem channel on slack. Register here: https://slackinvite.julialang.org/

Thanks.

[natschil]

@KristofferC: Thanks! The fact that I wanted to compute global coordinates was actually why I was looking at the reinit! function. I'll have a look at the reference, and if I still have further questions I'll ask them on slack.

[mohamed82008]

That slack invitation is really taking its time, is it ok to ask here until I can ask there?

[KristofferC]

I thought it would automatically accept any applicant at this point. Sure, go ahead :)

[mohamed82008]

What is the purpose of sorting during assembly?

[KristofferC]

To be able to go through the column in the global stiffness matrix linearly and just "tick off" one dof after the other. I found this usually faster than doing the binary search for each dof.

[mohamed82008]

Hm, slick! I guess it is faster because of the maxlookups.

[mohamed82008]

I would also appreciate an answer to https://discourse.julialang.org/t/propagate-inbounds-and-boundscheck-vs-checkbounds/7147 thank you :)