ferrite-fem / ferrite.jl Goto Github PK

So that it is simple to get mesh files written for Matlab inte JuAFEM.

Package to use: https://github.com/simonster/MAT.jl

Some of the function names feels a bit awkward. I suggest that we change some of them so that they are easier to "say".

function_scalar_value -> scalar_function_value # "value of a scalar function" rather than "value of a function scalar"
function_vector_value -> vector_function_value
function_scalar_gradient -> scalar_function_gradient
...

Get nodes from elementset / boundaryset.

What is the best approach to assemble the matrix that arises from the following weak form on [-1,1]:

\int 2x v_x u - (1-x^2) v_x u_x dx

v is the trial function and u is the test function. Subscripts are derivatives wrt x
I can set up the v_x u term and the v_x u_x term by using shape_value and shape_gradient, but need some guidance for how to include the 2x and 1-x^2 terms.

My assembler function so far:

function doassemble{dim}(cellvalues::CellScalarValues{dim}, K::SparseMatrixCSC, dh::DofHandler)
    # Trying the set up matrix for the linear form
    # \int u \mathcal L v\,dx
    # = \int 2x(b-a)∂v(x)*u(x) - b(1-x^2)∂v(x)∂u(x) dx

    assembler = start_assemble(K)

    n_basefuncs = getnbasefunctions(cellvalues)
    global_dofs = zeros(Int, ndofs_per_cell(dh))
    Ke = zeros(n_basefuncs, n_basefuncs) # Local FEM matrix

    @inbounds for (cellcount, cell) in enumerate(CellIterator(dh))
        fill!(Ke, 0)

        reinit!(cellvalues, cell)

        for q_point in 1:getnquadpoints(cellvalues)
            dΩ = getdetJdV(cellvalues, q_point)
            for i in 1:n_basefuncs
                vi = shape_value(cellvalues, q_point, i)
                dvi = shape_gradient(cellvalues, q_point, i)[1]
                for j in 1:n_basefuncs
                    vj = shape_value(cellvalues, q_point, j)
                    dvj = shape_gradient(cellvalues, q_point, j)[1]                   
                    Ke[i,j] += (dvi * vj) * dΩ # Missing 2x(b-a)
                    Ke[i,j] -= (dvi * dvj) * dΩ # Missing b(1-x^2)
                end
            end
        end

        celldofs!(global_dofs, cell)
        assemble!(assembler, global_dofs, Ke)
    end
end

Explain how the elements functions are being generated.

https://github.com/KristofferC/JuAFEM.jl/blob/7f5ae3f1e614f6213c8e7b8216d43a96172739c6/src/Grid/grid.jl#L124

I believe the way the all keyword is implemented is not correct here. This is because if all is false, then pass is false and there is no way for it to become true, so no cells will be added. Initializing pass to true would probably fix this, PR in the making.

This should just be handled by sending in a closure f(x) = f(x, t, ...).

Probably due to @nloops and friends? Mostly on v0.6, but pretty slow on v0.5 aswell:

v0.5:

julia> using JuAFEM

julia> @time qr = QuadratureRule{2, RefCube}(:legendre, 2);
  3.661755 seconds (5.95 M allocations: 239.781 MB, 2.25% gc time)

v0.6:

julia> using JuAFEM

julia> @time qr = QuadratureRule{2, RefCube}(:legendre, 2);
 10.056088 seconds (11.82 M allocations: 544.792 MiB, 1.60% gc time)

In the following setting

dim = 2
grid = generate_grid(Quadrilateral, (1, 1))
ip = Lagrange{dim, RefCube, 1}()
qr = QuadratureRule{dim, RefCube}(2)
cellvalues = CellScalarValues(qr, ip)

results in

JuAFEM.CellScalarValues{2,Float64,JuAFEM.RefCube} with 9 shape functions and 4 quadrature points

so i would expect the DofHandler to produce a 9x9 sparsity pattern, but it doesn't have the information for recognizing that.

https://github.com/JuliaCI/PkgBenchmark.jl

Document everything nicely.

Somehow, the Julia 0.6 package manager doesn't seem to get that the latest master is made for Julia 0.7/1.0 only. Is it possible to have a tagged version that works with Julia 0.6? The complaint is that Julia doesn't know the library LinearAlgebra.

I don't think it works .

Once the JuAFEM upgrade to Julia 0.7 is through, how would one include it in a non-registered package that uses REQUIRE? Just mentioning JuAFEM there seems to give an error.

It looks like the following vars are exported without defined references:

• CellFace
• FaceIterator

see:

https://github.com/KristofferC/JuAFEM.jl/blob/70447343a5b74360f330f4581a1fdfc3df22e0b3/src/exports.jl#L54

https://github.com/KristofferC/JuAFEM.jl/blob/70447343a5b74360f330f4581a1fdfc3df22e0b3/src/exports.jl#L96

Hi everyone,

I am a new learner of Julia and JuAFEM. I want to ask a preliminary question about calling the JuAFEM functions: In the Julia-REPL, I run some functions like that

julia> using JuAFEM, SparseArrays

julia> grid = generate_grid(Quadrilateral, (20, 20));

julia> dim = 2;

julia> ip = Lagrange{dim, RefCube, 1}();

I want to check the functions in interpolations.jl, first I run the getnbasefunctions() function successfully

julia> getnbasefunctions(ip)
4

but the getdim() function failed

julia> getdim(ip)
ERROR: UndefVarError: getdim not defined
Stacktrace:
 [1] top-level scope at none:0

the calling of value() function is also failed.

Actually, getnbasefunctions(), getdim(), value() are in the same file interpolations.jl, and getnbasefunctions() is exported by the exports.jl.

But in the cell_values.jl I find that

@assert getdim(func_interpol) == getdim(geom_interpol)
...
ξ -> value(geom_interpol, basefunc, ξ)

I'm wondering to know that why the functions in cell_values.jl can call the functions in interpolations.jl which are not exported by the exports.jl?

The function doassemble in the heat_equation.jl example has many julia comments in the code, but uses # rather than ##, so the code gets split up in the markdown presentation:
http://kristofferc.github.io/JuAFEM.jl/latest/examples/generated/heat_equation.html

I'm just checking whether that's the intention, or due to the change in how comments are processed in Literate.jl

The code for reinit! for FEBoundaryValues is very similar to the one for FECellValues. Just need to set the active boundary and then the rest is the same so should be able to pull it out to a common function.

Or something else. Ideas @fredrikekre?

For example, you might want a new thread to create a new assembler. So you call start_assemble which will reset the values in the sparse matrix which is not what you want.

https://github.com/KristofferC/JuAFEM.jl/pull/232/files is bugged because of this.

Only Float64 precision is fully supported. I will make a PR for this.

Done:

Elements

Springs

• spring1e
• spring1s

Bars

• bar2e
• bar2s
• bar2g

Solids (TODO: plane stress)

• plante
• plants
• plantf

• plani4e
• plani4s
• plani4f

• plani8e
• plani8s
• plani8f

• soli8e
• soli8s
• soli8f

Heat flow

• flw2te
• flw2ts

• flw2i8e
• flw2i8s

• flw2i4e
• flw2i4s

• flw2i8e
• flw2i8s

Materials

• hooke

Utilities

• assem
• solveq
• extract
• statcon
• coordxtr

Graphics

• eldraw2
• eldisp2

To do:

Elements

Bars

• bar3e
• bar3s

Beams

• beam2e
• beam2s

• beam2t
• beam2ts

• beam2w
• beam2ws

• beam2g
• beam2gs

• beam2d
• beam2ds

• beam3e
• beam3s

Materials

• mises
• dmises

Utilities

• eigen
• insert

Grahpics

• fill
• eldia2
• elflux2
• elprinc2
• scalfact2
• pltscalb2

Would be good to have.

So that for a 2Dcase:

g = zeros(6) # residual
for gp in 1:nqp
    for i in 1:ndofs
        g[i] += (shape_symmetric_gradient(fev, qp, i) : σ) dΩ
    end
end

This line https://github.com/KristofferC/JuAFEM.jl/blob/a8d976d4e9d4bba415c3d2dbf062a75f8d814b2c/src/Dofs/DofHandler.jl#L184 will call a function that is not defined for faces with 4 corner vertices, and that have dofs on the inside of the face. This will be the case when using a third order Hexahedron with dofs on the inside of the face, pretty niche but technically a bug!

As far as I can tell the current generate_grid() for Tetrahedral elements in 3D results in discontinuous "hat" functions.

To see this, notice that for every hexahedron the line (1,0,0) to (1,1,1) passes through a face and separates two elements. But for the adjacent hexahedron, this line is within an element. See also http://www.baumanneduard.ch/Splitting%20a%20cube%20in%20tetrahedras2.htm (which is cited in generate_grid which does a rotated version of #13): if you look at the image of #13, you'll see that the lines on adjacent cubes don't fit together.

For a 3D case, the FEBoundaryValues should be a 2D space, embedded in a 3D space. Like dealii has a dim and a spacedim.

https://github.com/KristofferC/JuAFEM.jl/blob/a8d976d4e9d4bba415c3d2dbf062a75f8d814b2c/src/Dofs/ConstraintHandler.jl#L211

The function faces only returns tuples of corner points, so the boundary conditions on non-corner points will not be written.

After learning how to use/using JuAFEM for about 2 weeks, I must say that it is a great tool, very understandable and hackable design, and it seems general enough to be expanded upon greatly!

I've used JuAFEM basically as a discretization tool for calculating a continuum model's total energy and gradients/hessians, as discussed on juafem-forwarddiff-nlsolve-finite-elements-for-nonlinear-problems It performs very well, vastly outperforming other implementations that colleagues of mine have done in c++ or Mathematica (duh).

I'll now give my wishlist for features as requested in the thread:

• Easier ways to access the fields inside a DofHandler, and create starting conditions e.g.:

function startingconditions!(dofvec, dh)
    mid = div(ndofs_per_cell(dh), 2)
    for cell in CellIterator(dh)
        globaldofs = celldofs(cell)
        for i=1:3:mid
            it = div(i, 3) + 1
            dofvec[globaldofs[mid+i]] = -0.53
            dofvec[globaldofs[mid+i+1]] = 0.53tanh(cell.coords[it][1]/5)
            dofvec[globaldofs[mid+i+2]] = 0.53
        end
    end
end

In light of the other design, something like fielddofs(:Field, cell) similar to celldofs(cell), seems a good way. Also fielddofs(:Field, dh::DofHandler) would be very nice to calculate certain things from the finalized result.
Edit: see #187 for implementation

• A way to add fancier constraints, such as bounds on fields etc

• More options to generate grids/read them from formats, although I think this is already somewhat possible with various other packages. What would be cool is if it was possible to generate non uniform grids, defining some kind of rule for element size/mesh density, depending on the coordinates.

• Higher flexibility on combining field dimensions with geometry dimensions, e.g. 3D field on a 2D geometry.

• A way to thread over the CellIterator loop.
  This seems a bit tricky since assembling residus or stiffness tensors is not trivial multithreaded. It really depends on the grid what needs to be done. The way I do it now, which is very much not scalable is pretty brutal (I have two fields, :u, :P, and thus a dofvector uP):

mutable struct ModelCaches{T, CV, CI, CF}
    ∇fs::Vector{Vector{T}}
    cvus::Vector{CV}
    cvPs::Vector{CV}
    cellits::CI
    configs::CF
end

function ModelCaches(∇f::Vector{T}, cvu, cvP, dh) where T
    cellits = [CellIterator(dh) for t=1:nth]
    ellen   = length(celldofs(cellits[1]))
    ∇fs     = [zeros(∇f) for t=1:nth]
    cvus    = [copy(cvu) for t=1:nth]
    cvPs    = [copy(cvP) for t=1:nth]
    configs = [GradientConfig(nothing, zeros(T, ellen), Chunk(ellen)) for t=1:nth]
    return ModelCaches(∇fs, cvus, cvPs, cellits, configs)
end

function ∇F!(∇f::Vector{T} , uP::Vector{T}, caches, params) where T
    @threads for i=1:length(caches.∇fs)
        fill!(caches.∇fs[i], zero(T))
    end
    results = [zeros(T, length(celldofs(caches.cellits[1]))) for t=1:nth]
    @threads for i=1:length(caches.cellits[1])
        tid = threadid()
        reinit!(caches.cellits[tid], i)
        reinit!(caches.cvus[tid], caches.cellits[tid])
        reinit!(caches.cvPs[tid], caches.cellits[tid])
        ForwardDiff.gradient!(results[tid],
            uel -> element_potential(uel, caches.cvus[tid], caches.cvPs[tid], params),
            uP[celldofs(caches.cellits[tid])],
            caches.configs[tid])
        assemble!(caches.∇fs[tid], celldofs(caches.cellits[tid]), results[tid])
    end
    ∇f .= caches.∇fs[1]
    for fs in caches.∇fs[2:end]
        ∇f .+= fs
    end
end

This works, but is very much not scalable, and quite a lot of redundant things happen.

For now this is what I ran into the most. If I have more I'll update this issue. I could also maybe try to implement a couple of these things myself, no promises on that :p.
Anyway keep up the good work, and again thanks for helping me in the discourse thread!

After #159 we add boundary conditions on facesets. It is sometimes needed to lock individual nodes, to hinder e.g. rotation. We should have a way to do this.

It's annoying to have all the CALFEM stuff in the same repo as the other new stuff. I suggest moving all the CALFEM stuff to its own repo, CALFEM.jl perhaps, and just let that one use JuAFEM like it does now. This would make it easier to keep the examples using the new JuAFEM stuff and the CALFEM stuff separate.

Not sure if this is a bug in julia or a syntax error in the package. Anyway the following code confuses julia steprange and colon from ternary operator.

Code to reproduce the error:

using JuAFEM
grid = generate_grid(Tetrahedron, (20, 20, 20))
dh = DofHandler(grid)
ip = Lagrange{3, RefTetrahedron, 2}()

push!(dh, :u, 1, ip)
close!(dh)

line of the error:

https://github.com/KristofferC/JuAFEM.jl/blob/a007e3d0a650954472fe059150088abc69d567c6/src/Dofs/DofHandler.jl#L162

Code to fix error:
for edgedof in (dir == olddir ? (1:interpolation_info.nedgedofs) : (interpolation_info.nedgedofs:-1:1))

Cheers!

• Move to HTML
• Fixup name changes

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

Umbrella issue to track things to tweak:

• some way for a node to know it is on the boundary. Useful for something like this: addnodeset!(grid, "all nodes on boundary", onboundary)
  Edit: We can do this with dispatch on onboundary
• currently push! is used to add fields to a dofhandler, while for grid we have addnodeset/addcellset. Should we have addfield!() to be more consistent? (and for BC we simply have add!)
• ~ export some more stuff: free_dofs, dirichlet_dofs ? ~
• ~ reduce(K::Matrix, dbc) that reduces the matrix ~

There is no easy way right now to compute stuff needed for, e.g. surface traction. dealii has a FEFaceValues https://dealii.org/8.4.0/doxygen/deal.II/classFEFaceValues.html which has a reinit! method that takes a cell and a face number.

Not sure how this is done best.

Now that we will apply boundary conditions to facesets (#159) we need a addfaceset!(::Grid, ::String, ::Function) method similar to addnodeset! and addcellset!.

Since #71 this is not so clear IMO. Call the geometric interpolation functions for shape functions (since they are connected to the shape, e.g. geometry) and the function interpolation functions for base functions (they are the base for the FE-approximation)

This would probably mean that some current functions have to be renamed, e.g.

shape_value -> base_value
shape_gradient -> base_gradient
...

I am trying to write a model extractor using FreeCAD into JuAFEM and I ran into a problem. Some of the boundary conditions are defined over faces and the way faces are handled in JuAFEM is different from the way they are handled in FreeCAD. In FreeCAD, a face is just an element with a bunch of nodes, while in JuAFEM, a face is a local property of the cell with a local index in each cell. This means that a mechanism to match boundary faces is needed to read the boundary conditions data into JuAFEM from other files like .inp using FreeCAD.

Initially in https://github.com/KristofferC/JuAFEM.jl/pull/194, I had thought that for Dirichlet boundary conditions, it is possible to assign the BC on the nodes of the faces instead. But this is not correct since the field interpolation can have more support points than face nodes. So the only correct way to do this is to actually match the faces. Even more is the case for Neumann boundary conditions, where a traction force is applied on a face, and cannot be easily distributed among the nodes. Also ideally, the translation done back and forth should not result in any loss of information. So I have been thinking of the best way to make it easier to match faces in JuAFEM with faces in FreeCAD and I thought of the following.

If we cache the set of face ids that a boundary node is part of, in the form of (cell_id, local_face_id), then it is possible using the nodes' ids of the face on which the BC is applied on, to find the intersection of the sets and identify the matching face in JuAFEM. Caching a vector of boundary node ids, and another vector of sets in the grid struct is an idea, probably not the most cache efficient way though.

Please let me know what you think.

Is there a reason for throwing an error for negative Jacobian determinants instead of taking the absolute value? I've always seen the change of variables formula with an absolute value, but I could be missing something subtle where the orientation matters.

https://github.com/KristofferC/JuAFEM.jl/blob/efa2eff1eef971eeb077ef12cd4933c02a62a2a3/src/FEValues/cell_values.jl#L156-L157

https://github.com/stevengj/NBInclude.jl

For example:

function_space = Lagrange{2, RefTetrahedron, 1}()
quad_rule = GaussQuadrature(Dim{2}, RefTetrahedron(), 1)

Better would probably be

function_space = Lagrange(Dim{2}, RefTetrahedron(), 1)

so they look the same.

Edit: It doesn't work because then we need like Order{1}...

This is a petition to remove the unnecessary first row of the edof matrix. It is never used and makes life harder than it needs to be.

Advantages:

• Everywhere the edof matrix is used, the first row is removed
• Simplifies coding since, you can write edof(:,i) instead of edof(2:end,i)
• The code looks nicer, i.e edof(:,i) looks much better than edof(2:end,i)

Disadvantages:

• None (except some rewriting of code

Please vote here!

Nature is pleased with simplicity. And nature is no dummy
― Isaac New

I propose we change the QuadratureRule constructor to use dim and refshape as parameters.

Now:

julia> QuadratureRule(Dim{2}, RefCube(), 2)
JuAFEM.QuadratureRule{2,Float64}([1.0,1.0,1.0,1.0],ContMechTensors.Tensor{1,2,Float64,2}[[-0.57735,-0.57735],[0.57735,-0.57735],[-0.57735,0.57735],[0.57735,0.57735]])

Desired:

julia> QuadratureRule{2, RefCube}(2)
JuAFEM.QuadratureRule{2,Float64}([1.0,1.0,1.0,1.0],ContMechTensors.Tensor{1,2,Float64,2}[[-0.57735,-0.57735],[0.57735,-0.57735],[-0.57735,0.57735],[0.57735,0.57735]])

will have conflicts with the changes in #73 I think so maybe change after that is merged. Or I might include it in that PR, if it helps with the coding of FEFaceValues :)