Skip to content

Commit

Permalink
Merge 3601116 into 86dd068
Browse files Browse the repository at this point in the history
  • Loading branch information
@KristofferC
KristofferC committed Dec 10, 2018
2 parents 86dd068 + 3601116 commit aff98b3
Show file tree
Hide file tree
Showing 25 changed files with 310 additions and 378 deletions.
6 changes: 6 additions & 0 deletions .travis.yml
View file
Expand Up @@ -8,5 +8,11 @@ julia:
- 1.0
- nightly

script:
- julia -e 'using Pkg; Pkg.add(PackageSpec(name="Tensors", rev="master"));
Pkg.develop(PackageSpec(path=pwd()));
Pkg.build("FEMBasis");
Pkg.test("FEMBasis"; coverage=true)'

after_success:
- julia -e 'cd(Pkg.dir("FEMBasis")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())'
88 changes: 35 additions & 53 deletions README.md
View file
Expand Up @@ -13,10 +13,10 @@ code. As a concrete example, to generate basis functions for a standard 10-node
tetrahedron one can write

```julia
code = create_basis(
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
Expand All @@ -27,30 +27,30 @@ code = create_basis(
(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
```julia
begin
mutable struct Tet10
end
function Base.size(::Type{Tet10})
return (3, 10)
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
function Base.length(::Type{Tet10})
return 10
end
Base.@pure function Base.length(::Type{Tet10})
return 10
end
function FEMBasis.get_reference_element_coordinates(::Type{Tet10})
return ((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))
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::Matrix{T}, xi::Tuple{T, T, T}) where T
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
Expand All @@ -66,39 +66,20 @@ begin
end
return N
end
function FEMBasis.eval_dbasis!(::Type{Tet10}, dN::Matrix{T}, xi::Tuple{T, T, T}) where T
function FEMBasis.eval_dbasis!(::Type{Tet10}, dN::Vector{<:Vec{3}}, xi::Vec)
@assert length(dN) == 10
(u, v, w) = xi
begin
dN[1, 1] = -3.0 + 4.0v + 4.0w + 2.0 * (2u)
dN[2, 1] = -3.0 + 4.0u + 4.0w + 2.0 * (2v)
dN[3, 1] = -3.0 + 4.0v + 4.0u + 2.0 * (2w)
dN[1, 2] = -1 + 2.0 * (2u)
dN[2, 2] = 0
dN[3, 2] = 0
dN[1, 3] = 0
dN[2, 3] = -1 + 2.0 * (2v)
dN[3, 3] = 0
dN[1, 4] = 0
dN[2, 4] = 0
dN[3, 4] = -1 + 2.0 * (2w)
dN[1, 5] = 4.0 + -4.0v + -4.0w + -4.0 * (2u)
dN[2, 5] = -4.0u
dN[3, 5] = -4.0u
dN[1, 6] = 4.0v
dN[2, 6] = 4.0u
dN[3, 6] = 0
dN[1, 7] = -4.0v
dN[2, 7] = 4.0 + -4.0u + -4.0w + -4.0 * (2v)
dN[3, 7] = -4.0v
dN[1, 8] = -4.0w
dN[2, 8] = -4.0w
dN[3, 8] = 4.0 + -4.0v + -4.0u + -4.0 * (2w)
dN[1, 9] = 4.0w
dN[2, 9] = 0
dN[3, 9] = 4.0u
dN[1, 10] = 0
dN[2, 10] = 4.0w
dN[3, 10] = 4.0v
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
Expand All @@ -109,23 +90,23 @@ Also more unusual elements can be defined. For example, pyramid element cannot b
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:
```julia
code = create_basis(
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)
```
Expand All @@ -138,15 +119,16 @@ 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:

```julia
using Tensors
B = Quad4()
X = ([0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0])
u = ([0.0, 0.0], [1.0, -1.0], [2.0, 3.0], [0.0, 0.0])
grad(B, u, X, (0.0, 0.0))
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
```julia
2×2 Array{Float64,2}:
2×2 Tensors.Tensor{2,2,Float64,4}:
1.5 0.5
1.0 2.0
```
Expand Down
1 change: 1 addition & 0 deletions REQUIRE
View file
@@ -1,2 +1,3 @@
julia 0.7
Calculus
Tensors
16 changes: 14 additions & 2 deletions src/FEMBasis.jl
View file
Expand Up @@ -3,9 +3,22 @@

module FEMBasis

import Calculus
using LinearAlgebra
using Tensors
# Reexport Vec
export Vec

abstract type AbstractBasis end
const Vecish{N, T} = Union{NTuple{N, T}, Vec{N, T}}

abstract type AbstractBasis{dim} end
# Forward methods on instances to types
Base.length(B::T) where T<:AbstractBasis = length(T)
Base.size(B::T) where T<:AbstractBasis = size(T)
eval_basis!(B::T, N, xi) where T<:AbstractBasis = eval_basis!(T, N, xi)
eval_dbasis!(B::T, dN, xi) where T<:AbstractBasis = eval_dbasis!(T, dN, xi)
eval_basis(B::AbstractBasis{dim}, xi) where {dim} = eval_basis!(B, zeros(length(B)), xi)
eval_dbasis(B::AbstractBasis{dim}, xi) where {dim} = eval_dbasis!(B, zeros(Vec{dim}, length(B)), xi)

include("subs.jl")
include("vandermonde.jl")
Expand Down Expand Up @@ -35,7 +48,6 @@ include("nurbs_surface.jl")
export NSurf
include("nurbs_solid.jl")
export NSolid

include("math.jl")
export interpolate, interpolate!, jacobian
export grad, grad!
Expand Down
55 changes: 27 additions & 28 deletions src/create_basis.jl
View file
@@ -1,10 +1,6 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE

using Calculus

import Base.parse

function get_reference_element_coordinates end
function eval_basis! end
function eval_dbasis! end
Expand All @@ -13,16 +9,18 @@ function calculate_interpolation_polynomials(p, V)
basis = []
first(p.args) == :+ || error("Use only summation between terms of polynomial")
args = p.args[2:end]
N = size(V, 1)
b = zeros(N)
for i in 1:N
n = size(V, 1)
b = zeros(n)
for i in 1:n
fill!(b, 0.0)
b[i] = 1.0
# TODO: Think about numerical stability with
# inverting Vandermonde matrix?
solution = V \ b
N = Expr(:call, :+)
for (ai, bi) in zip(solution, args)
isapprox(ai, 0.0) && continue
push!(N.args, simplify( :( $ai * $bi ) ))
push!(N.args, Calculus.simplify( :( $ai * $bi ) ))
end
push!(basis, N)
end
Expand All @@ -31,32 +29,32 @@ end

function calculate_interpolation_polynomial_derivatives(basis, D)
vars = [:u, :v, :w]
dbasis = []
for N in basis
dbasis = Matrix(undef, D, length(basis))
for (i, N) in enumerate(basis)
partial_derivatives = []
for j in 1:D
push!(partial_derivatives, simplify(differentiate(N, vars[j])))
dbasis[j, i] = Calculus.simplify(Calculus.differentiate(N, vars[j]))
end
push!(dbasis, partial_derivatives)
end
return dbasis
end

function create_basis(name, description, X::NTuple{N, NTuple{D, T}}, p::Expr) where {N, D, T}
function create_basis(name, description, X::Vector{<:Vecish{D}}, p::Expr) where D
@debug "create basis given antsatz polynomial" name description X p
V = vandermonde_matrix(p, X)
basis = calculate_interpolation_polynomials(p, V)
return create_basis(name, description, X, basis)
end

function create_basis(name, description, X::NTuple{N, NTuple{D, T}}, basis) where {N, D, T}
function create_basis(name, description, X::Vector{<:Vecish{D}}, basis::Vector) where D
@assert length(X) == length(basis)
@debug "create basis given basis functions" name description X basis
dbasis = calculate_interpolation_polynomial_derivatives(basis, D)
return create_basis(name, description, X, basis, dbasis)
return create_basis(name, description, Vec.(X), basis, dbasis)
end

function create_basis(name, description, X::NTuple{N, NTuple{D, T}}, basis, dbasis) where {N, D, T}

function create_basis(name, description, X::Vector{<:Vecish{D, T}}, basis, dbasis) where {D, T}
N = length(X)
@debug "create basis given basis functions and derivatives" name description X basis dbasis

Q = Expr(:block)
Expand All @@ -66,9 +64,7 @@ function create_basis(name, description, X::NTuple{N, NTuple{D, T}}, basis, dbas

V = Expr(:block)
for i=1:N
for k=1:D
push!(V.args, :(dN[$k,$i] = $(dbasis[i][k])))
end
push!(V.args, :(dN[$i] = Vec(float.(tuple($(dbasis[:, i]...))))))
end

if D == 1
Expand All @@ -80,10 +76,10 @@ function create_basis(name, description, X::NTuple{N, NTuple{D, T}}, basis, dbas
end

code = quote
struct $name <: FEMBasis.AbstractBasis
struct $name <: FEMBasis.AbstractBasis{$D}
end

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

Expand All @@ -92,27 +88,30 @@ function create_basis(name, description, X::NTuple{N, NTuple{D, T}}, basis, dbas
j == 2 && return $N
end

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

function FEMBasis.get_reference_element_coordinates(::Type{$name})
return $X
end

function FEMBasis.eval_basis!(::Type{$name}, N::Matrix{T}, xi) where T
@inline function FEMBasis.eval_basis!(::Type{$name}, N::Vector{<:Number}, xi::Vec)
@assert length(N) == $N
$unpack
$Q
@inbounds $Q
return N
end

function FEMBasis.eval_dbasis!(::Type{$name}, dN::Matrix{T}, xi) where T
@inline function FEMBasis.eval_dbasis!(::Type{$name}, dN::Vector{<:Vec{$D}}, xi::Vec)
@assert length(dN) == $N
$unpack
$V
@inbounds $V
return dN
end
end

return code
end

create_basis_and_eval(args...) = eval(create_basis(args...))

0 comments on commit aff98b3

Please sign in to comment.