State of JuAFEM

Specific Domains Numerics
, ,
1

I’m reading the JuAFEM dev docs and it’s not immediately apparent whether a few things are supported. Are:

  • general unstructured grids (tet, hex)
  • wedge elements
  • nonconforming and adaptive meshes

supported?

Secondly, it seems that I could use the cell iterators and face iterators to implement discontinuous Galerkin schemes – is that correct? Or do the assembly routines assume continuous Galerkin linear systems?

I’m considering implementing some DG and HDG schemes in the library and it looks super promising since I wouldn’t have to code all the basis functions and master element information from scratch!

2

general unstructured grids (tet, hex)

Yes, there is no assumption about any structure in the grid.

wedge elements

There are no built-in “interpolations” for wedge elements but it should be fairly easy to add (as an example here is one implementation of an interpolation: https://github.com/KristofferC/JuAFEM.jl/blob/20d14ac223598c926d13023a511a48f0c5c1ce33/src/interpolations.jl#L176-L208.)

  • nonconforming and adaptive meshes

Not really, the Grid type we have in JuAFEM is quite simple and there isn’t really any hanging node support implemented.

I would note that JuAFEM is (at least intended) to be more a library than a “framework”, that is, you should be able to pick out the parts of JuAFEM that you like but not be forced “to use everything of JuAFEM just to use some of JuAFEM”. Just to give some example of toying around in the REPL:

julia> dim = 2; element = RefCube; order = 2
2

julia> ip = JuAFEM.Lagrange{dim, element, order}()
Lagrange{2,RefCube,2}()

julia> JuAFEM.reference_coordinates(ip) # coordinates of the reference element
9-element Array{Tensor{1,2,Float64,2},1}:
 [-1.0, -1.0]
 [1.0, -1.0]
 [1.0, 1.0]
 [-1.0, 1.0]
 [0.0, -1.0]
 [1.0, 0.0]
 [0.0, 1.0]
 [-1.0, 0.0]
 [0.0, 0.0]

julia> [JuAFEM.value(ip, i, Vec(0.5, 0.3)) for i in 1:9] # local shape functions in xi = [0.5, 0.3]
9-element Array{Float64,1}:
  0.013125
 -0.039375
  0.073125
 -0.024375
 -0.07875
  0.34125
  0.14625
 -0.11375
  0.6825

julia> sum(ans) # shape functions sum to 1
1.0

julia> qr = QuadratureRule{dim, element}(order); # typically want to integrate using some quadrature

julia> getweights(qr)
4-element Array{Float64,1}:
 0.9999999999999996
 0.9999999999999996
 0.9999999999999996
 0.9999999999999996

julia> getpoints(qr)
4-element Array{Tensor{1,2,Float64,2},1}:
 [-0.5773502691896257, -0.5773502691896257]
 [0.5773502691896257, -0.5773502691896257]
 [-0.5773502691896257, 0.5773502691896257]
 [0.5773502691896257, 0.5773502691896257]


julia> fe_values = CellScalarValues(qr, ip); # object used to create "global" shape values / function values etc

julia> JuAFEM.reinit!(fe_values, 5 .* JuAFEM.reference_coordinates(ip)); # initialize the object using the coordinates of some element, here just a bigger square than the reference one

julia> [shape_value(fe_values, #=q_point=# 1, i) for i in 1:9] # global shape values in the first quadrature point
9-element Array{Float64,1}:
  0.20733615597604868
 -0.055555555555555546
  0.014886066246173483
 -0.055555555555555546
  0.30356120084098637
 -0.08133897861876414
 -0.08133897861876414
  0.30356120084098637
  0.44444444444444453

julia> sum(ans)
1.0

etc.

2 Likes
3

Ok, thanks for that @kristoffer.carlsson Some follow up questions:

  1. What is the relationship/difference between JuAFEM and JuliaFEM?
  2. Where in the JuliaFEM packages can I find index maps that connect the DOFs for assembly to a linear system?

Thanks in advance!