juliafem / fembasis.jl Goto Github PK

FEMBasis contains interpolation routines for finite element function spaces. Given ansatz and coordinates of domain, shape functions are calculated symbolically in a very general way to get efficient code. Shape functions can also be given directly and in that case partial derivatives are calculated automatically.

Home Page: https://juliafem.github.io/FEMBasis.jl/latest

License: MIT License

Julia 100.00%

fem interpolation

fembasis.jl's Introduction

FEMBasis.jl contains interpolation routines for standard finite element function spaces. Given ansatz and coordinates of domain, interpolation functions are calculated symbolically in a very general way to get efficient code. As a concrete example, to generate basis functions for a standard 10-node tetrahedron one can write

code = FEMBasis.create_basis(
    :Tet10,
    "10 node quadratic tetrahedral element",
    [
     (0.0, 0.0, 0.0), # N1
     (1.0, 0.0, 0.0), # N2
     (0.0, 1.0, 0.0), # N3
     (0.0, 0.0, 1.0), # N4
     (0.5, 0.0, 0.0), # N5
     (0.5, 0.5, 0.0), # N6
     (0.0, 0.5, 0.0), # N7
     (0.0, 0.0, 0.5), # N8
     (0.5, 0.0, 0.5), # N9
     (0.0, 0.5, 0.5), # N10
    ],
    :(1 + u + v + w + u*v + v*w + w*u + u^2 + v^2 + w^2),
   )

The resulting code is

    struct Tet10 <: FEMBasis.AbstractBasis{3}
    end
    Base.@pure function Base.size(::Type{Tet10})
            return (3, 10)
        end
    function Base.size(::Type{Tet10}, j::Int)
        j == 1 && return 3
        j == 2 && return 10
    end
    Base.@pure function Base.length(::Type{Tet10})
            return 10
        end
    function FEMBasis.get_reference_element_coordinates(::Type{Tet10})
        return Tensors.Tensor{1,3,Float64,3}[[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0], [0.5, 0.0, 0.0], [0.5, 0.5, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 0.5], [0.5, 0.0, 0.5], [0.0, 0.5, 0.5]]
    end
    function FEMBasis.eval_basis!(::Type{Tet10}, N::Vector{<:Number}, xi::Vec)
        assert length(N) == 10
        (u, v, w) = xi
        begin
            N[1] = 1.0 + -3.0u + -3.0v + -3.0w + 4.0 * (u * v) + 4.0 * (v * w) + 4.0 * (w * u) + 2.0 * u ^ 2 + 2.0 * v ^ 2 + 2.0 * w ^ 2
            N[2] = -u + 2.0 * u ^ 2
            N[3] = -v + 2.0 * v ^ 2
            N[4] = -w + 2.0 * w ^ 2
            N[5] = 4.0u + -4.0 * (u * v) + -4.0 * (w * u) + -4.0 * u ^ 2
            N[6] = +(4.0 * (u * v))
            N[7] = 4.0v + -4.0 * (u * v) + -4.0 * (v * w) + -4.0 * v ^ 2
            N[8] = 4.0w + -4.0 * (v * w) + -4.0 * (w * u) + -4.0 * w ^ 2
            N[9] = +(4.0 * (w * u))
            N[10] = +(4.0 * (v * w))
        end
        return N
    end
    function FEMBasis.eval_dbasis!(::Type{Tet10}, dN::Vector{<:Vec{3}}, xi::Vec)
        @assert length(dN) == 10
        (u, v, w) = xi
        begin
            dN[1] = Vec(-3.0 + 4.0v + 4.0w + 2.0 * (2u), -3.0 + 4.0u + 4.0w + 2.0 * (2v), -3.0 + 4.0v + 4.0u + 2.0 * (2w))
            dN[2] = Vec(-1 + 2.0 * (2u), 0, 0)
            dN[3] = Vec(0, -1 + 2.0 * (2v), 0)
            dN[4] = Vec(0, 0, -1 + 2.0 * (2w))
            dN[5] = Vec(4.0 + -4.0v + -4.0w + -4.0 * (2u), -4.0u, -4.0u)
            dN[6] = Vec(4.0v, 4.0u, 0)
            dN[7] = Vec(-4.0v, 4.0 + -4.0u + -4.0w + -4.0 * (2v), -4.0v)
            dN[8] = Vec(-4.0w, -4.0w, 4.0 + -4.0v + -4.0u + -4.0 * (2w))
            dN[9] = Vec(4.0w, 0, 4.0u)
            dN[10] = Vec(0, 4.0w, 4.0v)
        end
        return dN
    end
end

Also more unusual elements can be defined. For example, pyramid element cannot be descibed with ansatz, but it's still possible to implement by defining shape functions, Calculus.jl is taking care of defining partial derivatives of function:

code = FEMBasis.create_basis(
    :Pyr5,
    "5 node linear pyramid element",
    [
     (-1.0, -1.0, -1.0), # N1
     ( 1.0, -1.0, -1.0), # N2
     ( 1.0,  1.0, -1.0), # N3
     (-1.0,  1.0, -1.0), # N4
     ( 0.0,  0.0,  1.0), # N5
    ],
    [
     :(1/8 * (1-u) * (1-v) * (1-w)),
     :(1/8 * (1+u) * (1-v) * (1-w)),
     :(1/8 * (1+u) * (1+v) * (1-w)),
     :(1/8 * (1-u) * (1+v) * (1-w)),
     :(1/2 * (1+w)),
    ],
   )
eval(code)

Basis function can have internal variables if needed, e.g. variable dof basis like hierarchical basis functions or NURBS.

It's also possible to do some very common FEM calculations, like calculate Jacobian or gradient of some variable with respect to some coordinates. For example, to calculate displacement gradient du/dX in unit square [0,1]^2, one could write:

using Tensors
B = Quad4()
X = Vec.([(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)])
u = Vec.([(0.0, 0.0), (1.0, -1.0), (2.0, 3.0), (0.0, 0.0)])
grad(B, u, X, Vec(0.0, 0.0))

Result is

2×2 Tensors.Tensor{2,2,Float64,4}:
 1.5  0.5
 1.0  2.0

fembasis.jl's People

Contributors

Stargazers

Watchers

Forkers

tkelman sudhakaryogaraj mlubin batmanabcdefg xtalax juliatagbot bocc haoxiangmiao hsugawa8651 mfkiwl standardgalactic

fembasis.jl's Issues

Documentation index is not descriptive

Documentation is empty. Now it looks like this.

Use expressions instead of strings?

Passing Julia expressions as strings into create_basis seems really odd. Is there a reason it's not just an expression?

Perhaps stop punning on Base functions?

FEMBasis does some "puns" on Base functions, for example:

        function Base.size(::Type{$name})
            return ($D, $N)
        end

        function Base.size(::Type{$name}, j::Int)
            j == 1 && return $D
            j == 2 && return $N
        end

        function Base.length(::Type{$name})
            return $N
        end

length is typically defined for collections (or at least when the instance of the type is a collection). AbstractBasis is not iterable and it probably shouldn't be. Instead of defining these functions I think they should have their own names ndim, nbasis.
Overloading Base functions are typically meant so you can write generic code but there is no way to write generic code with puns like this.

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Implement hierarchical basis elements

Implement P elements

API changes to FEMBasis 0.3

In the next minor, we should somehow solve the issue of having multiple basis functions for the same topology. Motivation comes here: https://github.com/JuliaFEM/FEMPlates.jl/blob/master/src/dkt.jl#L34

It looks that this is kind of normal situation, that some basis functions need to be calculated "in the fly" based on the actual geometry. And the same topology can have different basis functions depending on the problem. We can default the basis function evaluation to Lagrange nodal basis, but we should also offer an option to create more complex basis functions which cannot necessarily be constructed beforehand.

For constructing nodal basis functions, we could have something like

struct Tri3 end
struct Lagrange end
X = ((0.0,0.0), (1.0,0.0), (0.0,1.0))
p = :(1 + u + v)
@create_basis(Lagrange, Tri3, p, X)

And after that, we can evaluate basis functions

get_basis(Lagrange(), Tri3(), xi)
get_basis(Tri3(), xi) # defaults to Lagrange..?

Then, we could also implement e.g. that DKT basis, but for that, we need the geometry of the element:

struct DKT end
function FEMBasis.get_basis(::DKT, ::Tri3, X)
    # implementation
end
get_basis(DKT(), Tri3(), X, xi)

There is even more complicated situations, where the basis functions depends not only the element geometry, but also geometry from surrounding elements, like bi-orthogonal basis in contact mechanics.

Number of Warning messages in FEMBasis

FEMBasis seems to fill stdout log with warnings, f.ex.:
https://travis-ci.org/JuliaFEM/JuliaFEM.jl/jobs/428170944

Any idea where this roots from?

Fix coverage up to 100 %

https://coveralls.io/builds/15280039/source?filename=src%2Fcreate_basis.jl

Looks that Coverage.jl is not understanding the content inside macro.

Release version 0.3.0

This version uses Tensors.jl to get better performance.

@JuliaRegistrator register()

Release 0.3.2

Would it be possible to release version 0.3.2 to the General registry? Once this update is complete, several pieces of the ecosystem may begin working; c.f.:

JuliaFEM/JuliaFEM.jl#257
JuliaFEM/FEMBase.jl#70

Use StaticArrays.jl

It might be useful to use non-allocating static arrays all over the place

Implement Hermite etc. basis

We need to have methods to calculate basis with continuous derivatives. Resulting system of equations involves derivatives of functions. Concrete example p(x) = a0*1 + a1*x + a2*x^2 + a3*x^3, and equations are p(-1), p'(-1), p(1), p'(1). Calculus is able to automatically differentiate this.

"::Vec" Type declaration outdated?

Hello!

in lines 101 and 108 of create_basis.jl:

@inline function FEMBasis.eval_basis!(::Type{$name}, N::Vector{<:Number}, xi::Vec)
@inline function FEMBasis.eval_dbasis!(::Type{$name}, dN::Vector{<:Vec{$D}}, xi::Vec)

the type declaration throws me errors, when trying to use a Vector for xi. Is this on purpose or is it simply a relict from an older version? Changing them to xi::Vector fixed the issue. Furthermore, in the documentation of JuliaFEM.jl this is done with xi as a tuple.

This code should reproduce the example, with Julia v.1.9.2 and FEMBasis 0.3.2
B = Quad4()
N = zeros(1, length(B))
xi = [0.0, 0.0] #(or using (0.0, 0.0) as in JuliaFEM doc)
eval_basis!(B, N, xi)

I am trying to use FEMBasis as if i were using it from JuliaFEM, i.e. with the help of JuliaFEM documentation. A related issue is that basis functions are not defined in JuliaFEM, does anybody know about this?
using JuliaFEM
Tri3() #not defined error

I don't know if it is useful to open a seperate issue, or if this unavailability of the basis functions are related to this package. They seem to be exported in JuliaFEM, but not defined.

Thank you for the help in advance!

WARNING: eval from module __SANDBOX__ to FEMBasis:
:u = -1
  ** incremental compilation may be broken for this module **

WARNING: eval from module __SANDBOX__ to FEMBasis:
:v = -1
  ** incremental compilation may be broken for this module **

WARNING: eval from module __SANDBOX__ to FEMBasis:
:w = -1
  ** incremental compilation may be broken for this module **

WARNING: eval from module __SANDBOX__ to FEMBasis:
Expr(:call, :*, Expr(:call, :^, :u, 2)::Any, Expr(:call, :^, :v, 2)::Any, Expr(:call, :^, :w, 2)::Any)::Any
  ** incremental compilation may be broken for this module **

WARNING: eval from module __SANDBOX__ to FEMBasis:
:u = 1
  ** incremental compilation may be broken for this module **

WARNING: eval from module __SANDBOX__ to FEMBasis:
:v = -1
  ** incremental compilation may be broken for this module **

WARNING: eval from module __SANDBOX__ to FEMBasis:
:w = -1
  ** incremental compilation may be broken for this module **

WARNING: eval from module __SANDBOX__ to FEMBasis:
Expr(:call, :*, Expr(:call, :^, :u, 2)::Any, Expr(:call, :^, :v, 2)::Any, Expr(:call, :^, :w, 2)::Any)::Any
  ** incremental compilation may be broken for this module **

Examples are missing

We should generate some examples from test files using Literate.jl.

Finish NURBS elements

Calculate partial derivatives for elements