mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-20 01:59:59 +00:00
added normal calculation code to base
This commit is contained in:
Jukka Aho
File diff suppressed because one or more lines are too long
+223
-136
@@ -27,7 +27,7 @@ Several functions are inherited from Element abstract type:
|
||||
- get_number_of_basis_functions*
|
||||
- get_element_dimension *
|
||||
- get_basis *
|
||||
- get_dbasisdxi*
|
||||
- get_dbasisdxi *
|
||||
- get_dbasisdX
|
||||
- get_field
|
||||
- set_field
|
||||
@@ -36,46 +36,15 @@ Several functions are inherited from Element abstract type:
|
||||
|
||||
Which should work if element is defined following some rules. Functions marked with asterisk * are the ones which must necessarily to implement by your own.
|
||||
|
||||
=#
|
||||
|
||||
# These must be implemented for your own element
|
||||
get_number_of_basis_functions(el::Type{Element}) = nothing
|
||||
get_number_of_basis_functions(el::Element) = nothing
|
||||
get_element_dimension(el::Element) = nothing
|
||||
get_basis(el::Element, xi) = nothing
|
||||
get_dbasisdxi(el::Element, xi) = nothing
|
||||
get_connectivity(el::Element) = el.connectivity
|
||||
|
||||
"""
|
||||
Create new element with element_name to family element_family
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> @create_element(Seg2, CG, "2 node linear segment")
|
||||
"""
|
||||
macro create_element(element_name, element_family, element_description)
|
||||
# Logging.debug("Creating element ", element_name, ": ", element_description, "\n")
|
||||
eltype = esc(element_name)
|
||||
elfam = esc(element_family)
|
||||
quote
|
||||
global get_element_description
|
||||
type $eltype <: $elfam
|
||||
connectivity :: Array{Int, 1}
|
||||
fields :: Dict{Any, Any}
|
||||
end
|
||||
$eltype(connectivity) = $eltype(connectivity, Dict{Any, Any}())
|
||||
get_element_description(el::Type{$eltype}) = $element_description
|
||||
end
|
||||
end
|
||||
|
||||
#=
|
||||
|
||||
Start of example
|
||||
----------------
|
||||
|
||||
Example how to create new element. This is commented because I use code
|
||||
This is example how to create new element. This is commented because I use code
|
||||
generation for simple elements like Lagrage elements. Feel free to use
|
||||
code generation but elements can be of course created manually too!
|
||||
|
||||
abstract CG <: Element # create new element family "Continous Galerkin"
|
||||
|
||||
type Quad4 <: CG
|
||||
connectivity :: Array{Int, 1}
|
||||
fields :: Dict{Any, Any}
|
||||
@@ -110,96 +79,25 @@ End of example.
|
||||
|
||||
=#
|
||||
|
||||
# These must be implemented for your own element
|
||||
get_number_of_basis_functions(el::Type{Element}) = nothing
|
||||
get_number_of_basis_functions(el::Element) = nothing
|
||||
get_element_dimension(el::Element) = nothing
|
||||
get_basis(el::Element, xi) = nothing
|
||||
get_dbasisdxi(el::Element, xi) = nothing
|
||||
get_connectivity(el::Element) = el.connectivity
|
||||
|
||||
### LAGRANGE ELEMENTS ###
|
||||
include("lagrange.jl")
|
||||
|
||||
abstract CG <: Element # Lagrange (continous Galerkin) element family
|
||||
|
||||
"""
|
||||
Given polynomial P and coordinates of reference element, calculate
|
||||
Lagrange basis function and partial derivatives.
|
||||
"""
|
||||
function calculate_lagrange_basis(P, X)
|
||||
dim, nbasis = size(X)
|
||||
A = zeros(nbasis, nbasis)
|
||||
for i=1:nbasis
|
||||
A[i,:] = P(X[:, i])
|
||||
end
|
||||
# Logging.debug("Calculating inverse of A")
|
||||
invA = inv(A)'
|
||||
basis(xi) = invA*P(xi)
|
||||
dbasisdxi = ForwardDiff.jacobian(basis)
|
||||
basis, dbasisdxi
|
||||
end
|
||||
|
||||
"""
|
||||
Assign Lagrange basis for element.
|
||||
"""
|
||||
macro create_lagrange_basis(element_name, X, P)
|
||||
|
||||
# Logging.debug("Creating Lagrange basis for element ", element_name, ". ")
|
||||
eltype = esc(element_name)
|
||||
|
||||
quote
|
||||
|
||||
global get_number_of_basis_functions, get_element_dimension
|
||||
global get_basis, get_dbasisdxi
|
||||
|
||||
dim = size($X, 1)
|
||||
nbasis = size($X, 2)
|
||||
# Logging.debug("Number of basis functions: ", nbasis, ". ")
|
||||
# Logging.debug("Element dimension: ", dim)
|
||||
|
||||
get_number_of_basis_functions(el::Type{$(esc(element_name))}) = nbasis
|
||||
get_number_of_basis_functions(el::$(esc(element_name))) = nbasis
|
||||
get_element_dimension(el::$(esc(element_name))) = dim
|
||||
|
||||
basis, dbasisdxi = calculate_lagrange_basis($P, $X)
|
||||
get_basis(el::$eltype, xi) = basis(xi)
|
||||
get_dbasisdxi(el::$eltype, xi) = dbasisdxi(xi)
|
||||
# Logging.debug("Element ", $element_name, " created.")
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
# 0d Lagrange element
|
||||
|
||||
@create_element(Point1, CG, "1 node point element")
|
||||
|
||||
# 1d Lagrange elements
|
||||
|
||||
@create_element(Seg2, CG, "2 node linear line element")
|
||||
@create_lagrange_basis(Seg2, [-1.0 1.0], (xi) -> [1.0, xi[1]])
|
||||
|
||||
@create_element(Seg3, CG, "3 node quadratic line element")
|
||||
@create_lagrange_basis(Seg3, [-1.0 1.0 0.0], (xi) -> [1.0, xi[1], xi[1]^2])
|
||||
|
||||
# 2d Lagrange elements
|
||||
|
||||
@create_element(Quad4, CG, "4 node bilinear quadrangle element")
|
||||
@create_lagrange_basis(Quad4,
|
||||
[-1.0 1.0 1.0 -1.0
|
||||
-1.0 -1.0 1.0 1.0],
|
||||
(xi) -> [1.0, xi[1], xi[2], xi[1]*xi[2]])
|
||||
|
||||
# 3d Lagrange elements
|
||||
|
||||
@create_element(Tet10, CG, "10 node quadratic tetrahedron")
|
||||
@create_lagrange_basis(Tet10,
|
||||
[0.0 1.0 0.0 0.0 0.5 0.5 0.0 0.0 0.5 0.0
|
||||
0.0 0.0 1.0 0.0 0.0 0.5 0.5 0.0 0.0 0.5
|
||||
0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.5 0.5 0.5],
|
||||
(xi) -> [ 1.0, xi[1], xi[2], xi[3], xi[1]^2,
|
||||
xi[2]^2, xi[3]^2, xi[1]*xi[2], xi[2]*xi[3], xi[3]*xi[1]])
|
||||
|
||||
|
||||
### HIERARCHICAL ELEMENTS ###
|
||||
|
||||
### HIERARCHICAL P-ELEMENTS ###
|
||||
include("hierarchical.jl")
|
||||
|
||||
# Common element routines
|
||||
### COMMON ELEMENT ROUTINES ###
|
||||
|
||||
"""
|
||||
Test routine for element.
|
||||
Test routine for element. If this passes, element interface is properly
|
||||
defined.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
@@ -208,7 +106,7 @@ eltype::Type{Element}
|
||||
|
||||
Raises
|
||||
------
|
||||
This uses FactCheck and throws exception if element is not passing.
|
||||
This uses FactCheck and throws exceptions if element is not passing all tests.
|
||||
"""
|
||||
function test_element(eltype)
|
||||
Logging.info("Testing element $eltype")
|
||||
@@ -264,38 +162,46 @@ function test_element(eltype)
|
||||
Logging.info("Element $eltype passed tests.")
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
Get jacobian of element evaluated at point ξ on element in reference configuration.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
el::Element
|
||||
xi::Vector
|
||||
geometry_field::Any, optional
|
||||
|
||||
Returns
|
||||
-------
|
||||
Vector or Matrix
|
||||
depending on element type
|
||||
|
||||
Notes
|
||||
-----
|
||||
This function assumes that element has field :geometry defined.
|
||||
Big "J" comes from reference (undeformed) configuration.
|
||||
"""
|
||||
#function get_Jacobian(el::Element, xi)
|
||||
# dbasisdxi = get_dbasisdxi(el, xi)
|
||||
# X = get_field(el, :geometry)
|
||||
# J = X*dbasisdxi
|
||||
# return J
|
||||
#end
|
||||
function get_Jacobian(el::Element, xi, geometry_field=:geometry)
|
||||
function get_Jacobian(el::Element, xi, geometry_field=:Geometry)
|
||||
dinterpolate(el, geometry_field, xi)
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
Get jacobian of element evaluated at point ξ on element in current configuration.
|
||||
|
||||
Notes
|
||||
-----
|
||||
This function assumes that element has fields :geometry and :displacement defined.
|
||||
Small "j" comes from current (deformed) configuration.
|
||||
"""
|
||||
function get_jacobian(el::Element, xi)
|
||||
function get_jacobian(el::Element, xi, geometry_field=:Geometry, displacement_field=:displacement)
|
||||
dbasisdxi = get_dbasisdxi(el, xi)
|
||||
X = get_field(el, :geometry)
|
||||
u = get_field(el, :displacement)
|
||||
X = get_field(el, geometry_field)
|
||||
u = get_field(el, displacement_field)
|
||||
j = (X+u)*dbasisdxi
|
||||
return j
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
Evaluate partial derivatives of basis, dbasis/dX
|
||||
"""
|
||||
@@ -305,6 +211,7 @@ function get_dbasisdX(el::Element, xi)
|
||||
dbasisdxi*inv(J)
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
Evaluate partial derivatives of basis, dbasis/dx
|
||||
"""
|
||||
@@ -314,32 +221,38 @@ function get_dbasisdx(el::Element, xi)
|
||||
dbasisdxi*inv(j)
|
||||
end
|
||||
|
||||
|
||||
""" Set field variable. """
|
||||
function set_field(el::Element, field_name, field_value)
|
||||
el.fields[field_name] = field_value
|
||||
end
|
||||
|
||||
|
||||
""" Create new empty field of some type. """
|
||||
function new_field(el::Element, field_name, field_type)
|
||||
function new_field!(el::Element, field_name, field_type)
|
||||
el.fields[field_name] = field_type[]
|
||||
end
|
||||
|
||||
|
||||
""" Push to existing field. """
|
||||
function push_field!(el::Element, field_name, field_value)
|
||||
push!(el.fields[field_name], field_value)
|
||||
end
|
||||
|
||||
|
||||
""" Get field variable. """
|
||||
function get_field(el::Element, field_name)
|
||||
el.fields[field_name]
|
||||
end
|
||||
|
||||
"""Evaluate some field in point ξ on element using basis functions.
|
||||
|
||||
"""
|
||||
Evaluate some field in point ξ on element using basis functions.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
el :: Element
|
||||
field :: Union{ASCIIString, Symbol}
|
||||
field :: Any
|
||||
xi :: Vector
|
||||
|
||||
Returns
|
||||
@@ -375,3 +288,177 @@ function dinterpolate(el::Element, field, xi::Vector)
|
||||
end
|
||||
return sum([fld[i]*dbasis[i,:] for i in 1:length(fld)])
|
||||
end
|
||||
|
||||
"""
|
||||
calculate "local" normals in elements, in a way that
|
||||
n = Nᵢnᵢ gives some reasonable results for ξ ∈ [-1, 1]
|
||||
"""
|
||||
function calculate_normals!(el::Element, field_name=:Normals)
|
||||
new_field!(el, field_name, Vector)
|
||||
for xi in Vector[[-1.0], [1.0]]
|
||||
t = dinterpolate(el, :Geometry, xi)
|
||||
n = [0 -1; 1 0]*t
|
||||
n /= norm(n)
|
||||
push_field!(el, field_name, n)
|
||||
end
|
||||
end
|
||||
|
||||
"""
|
||||
Alter normal field such that normals of adjacent elements are averaged.
|
||||
"""
|
||||
function average_normals!(elements, normal_field=:Normals)
|
||||
d = Dict()
|
||||
for el in elements
|
||||
c = get_connectivity(el)
|
||||
n = get_field(el, normal_field)
|
||||
for (ci, ni) in zip(c, n)
|
||||
d[ci] = haskey(d, ci) ? d[ci] + ni : ni
|
||||
end
|
||||
end
|
||||
for (ci, ni) in d
|
||||
d[ci] /= norm(d[ci])
|
||||
end
|
||||
for el in elements
|
||||
c = get_connectivity(el)
|
||||
new_normals = [d[ci] for ci in c]
|
||||
set_field(el, normal_field, new_normals)
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
# FIXME: These two needs integration -- maybe not in elements.jl ..?
|
||||
"""
|
||||
Fit field s.t. || ∫ (Nᵢ(ξ)αᵢ - f(el, ξ)) dS || -> min!
|
||||
|
||||
Parameters
|
||||
----------
|
||||
f::Function
|
||||
Needs to take (el::Element, xi::Vector) as argument
|
||||
fixed_coeffs::Int[]
|
||||
These coefficients are not changed during fitting -> constrained optimizatio
|
||||
"""
|
||||
function fit_field!(el::Element, field, f, fixed_coeffs=Int[])
|
||||
w = [
|
||||
128/225,
|
||||
(332+13*sqrt(70))/900,
|
||||
(332+13*sqrt(70))/900,
|
||||
(332-13*sqrt(70))/900,
|
||||
(332-13*sqrt(70))/900]
|
||||
xi = Vector[
|
||||
[0.0],
|
||||
[ 1/3*sqrt(5 - 2*sqrt(10/7))],
|
||||
[-1/3*sqrt(5 - 2*sqrt(10/7))],
|
||||
[ 1/3*sqrt(5 + 2*sqrt(10/7))],
|
||||
[-1/3*sqrt(5 + 2*sqrt(10/7))]]
|
||||
n = get_number_of_basis_functions(el)
|
||||
fld = get_field(el, field)
|
||||
nfld = length(fld[1])
|
||||
#Logging.debug("dim of field $field: $nfld")
|
||||
|
||||
M = zeros(n, n)
|
||||
b = zeros(n, nfld)
|
||||
for i=1:length(w)
|
||||
detJ = get_detJ(el, xi[i])
|
||||
N = get_basis(el, xi[i])
|
||||
M += w[i]*N*N'*detJ
|
||||
fi = f(el, xi[i])
|
||||
for j=1:nfld
|
||||
b[:, j] += w[i]*N*fi[j]*detJ
|
||||
end
|
||||
end
|
||||
|
||||
coeffs = zeros(n)
|
||||
for j=1:nfld
|
||||
for k=1:n
|
||||
coeffs[k] = fld[k][j]
|
||||
end
|
||||
if length(fixed_coeffs) != 0
|
||||
# constrained problem, some coefficients are fixed
|
||||
N = Int[] # rest of coeffs
|
||||
S = Int[] # fixed coeffs
|
||||
for i = 1:n
|
||||
if i in fixed_coeffs
|
||||
push!(S, i)
|
||||
else
|
||||
push!(N, i)
|
||||
end
|
||||
end
|
||||
lhs = M[N,N]
|
||||
rhs = b[N,j] - M[N,S]*coeffs[S]
|
||||
coeffs[N] = lhs \ rhs
|
||||
else
|
||||
coeffs[:] = M \ b[:,j]
|
||||
end
|
||||
for k=1:n
|
||||
fld[k][j] = coeffs[k]
|
||||
end
|
||||
end
|
||||
set_field(el, field, fld)
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
Fit field s.t. || ∫ ∂/∂ξ(∑Nᵢ(ξ)αᵢ)f(el, ξ) dS || -> min!
|
||||
"""
|
||||
function fit_derivative_field!(el::Element, field, f, fixed_coeffs=Int[])
|
||||
w = [
|
||||
128/225,
|
||||
(332+13*sqrt(70))/900,
|
||||
(332+13*sqrt(70))/900,
|
||||
(332-13*sqrt(70))/900,
|
||||
(332-13*sqrt(70))/900]
|
||||
xi = Vector[
|
||||
[0.0],
|
||||
[ 1/3*sqrt(5 - 2*sqrt(10/7))],
|
||||
[-1/3*sqrt(5 - 2*sqrt(10/7))],
|
||||
[ 1/3*sqrt(5 + 2*sqrt(10/7))],
|
||||
[-1/3*sqrt(5 + 2*sqrt(10/7))]]
|
||||
n = get_number_of_basis_functions(el)
|
||||
fld = get_field(el, field)
|
||||
nfld = length(fld[1])
|
||||
#Logging.debug("dim of field $field: $nfld")
|
||||
|
||||
M = zeros(n, n)
|
||||
b = zeros(n, nfld)
|
||||
for i=1:length(w)
|
||||
detJ = get_detJ(el, xi[i])
|
||||
dNdxi = get_dbasisdxi(el, xi[i])
|
||||
dNdX = dNdxi / detJ
|
||||
M += w[i]*dNdX*dNdX'*detJ
|
||||
fi = f(el, xi[i])
|
||||
for j=1:nfld
|
||||
b[:, j] += w[i]*dNdX*fi[j]*detJ
|
||||
end
|
||||
end
|
||||
|
||||
coeffs = zeros(n)
|
||||
for j=1:nfld
|
||||
for k=1:n
|
||||
coeffs[k] = fld[k][j]
|
||||
end
|
||||
if length(fixed_coeffs) != 0
|
||||
#Logging.info("constrained problem, some coefficients are fixed")
|
||||
N = Int[] # rest of coeffs
|
||||
S = Int[] # fixed coeffs
|
||||
for i = 1:n
|
||||
if i in fixed_coeffs
|
||||
push!(S, i)
|
||||
else
|
||||
push!(N, i)
|
||||
end
|
||||
end
|
||||
lhs = M[N,N]
|
||||
rhs = b[N,j] - M[N,S]*coeffs[S]
|
||||
coeffs[N] = lhs \ rhs
|
||||
else
|
||||
coeffs[:] = M \ b[:,j]
|
||||
end
|
||||
for k=1:n
|
||||
fld[k][j] = coeffs[k]
|
||||
end
|
||||
end
|
||||
set_field(el, field, fld)
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
Reference in New Issue
Block a user