Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-08-06 04:21:33 +00:00
Code Issues Packages Projects Releases Wiki Activity

normal tangential coordinate system

Browse Source
This commit is contained in:
Jukka Aho
2015-12-12 10:08:20 +02:00
parent d960ec5824
commit 9013aae07c
7 changed files with 162 additions and 42 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/links.md
+10 -6
View File
@@ -1,12 +1,13 @@
Discretization
==============
Links
=====
Discretization
--------------
http://code.activestate.com/recipes/579021-delaunay-triangulation/
Interpolation
=============
-------------
- http://www.cs.rpi.edu/~flaherje/pdf/fea4.pdf
- http://www.sd.ruhr-uni-bochum.de/downloads/Shape_funct.pdf
- http://what-when-how.com/the-finite-element-method/fem-for-3d-solids-finite-element-method-part-1/
@@ -18,10 +19,13 @@ Interpolation
- http://www.researchgate.net/publication/267082822_Unified_isoparametric_3D_Lagrange_finite_elements
Integration
===========
-----------
- http://arxiv.org/pdf/1411.1341.pdf
Hierarchial shape functions
---------------------------
- https://www.math.vt.edu/people/adjerids/research/papers/basis.pdf
- edited by Ari
Solvers
-------
- https://github.com/JuliaSparse/MultiFrontalCholesky.jl
src/assembly.jl
+3 -1
View File
@@ -61,7 +61,9 @@ function assemble(problem::AllProblems, time::Real, nchunks=10)
for (j, elrange) in enumerate(slices)
sub_assembly = assemble(problem, elrange, time)
append!(assembly, sub_assembly)
info("Assembly: ", round(j/nchunks*100,1), " % done. ")
if ne > 100
info("Assembly: ", round(j/nchunks*100,1), " % done. ")
end
end
# optimize!(assembly)
# dim = length(assembly.stiffness_matrix.I)
src/elements.jl
+23
View File
@@ -210,3 +210,26 @@ function Base.haskey(element::Element, what)
haskey(element.fields, what)
end
""" Calculate local normal-tangential coordinates for element. """
function calculate_normal_tangential_coordinates!{E}(element::Element{E}, time::Real)
proj(u, v) = dot(v, u) / dot(u, u) * u
ntcoords = Matrix[]
refcoords = get_reference_element_coordinates(E)
x = element("geometry", time)
for xi in refcoords
dN = get_dbasis(E, xi)*x
normal = cross(dN[:,1], dN[:,2])
normal /= norm(normal)
u1 = normal
j = indmax(abs(u1))
v2 = zeros(3)
v2[mod(j,3)+1] = 1.0
u2 = v2 - proj(u1, v2)
u3 = cross(u1, u2)
tangent1 = u2/norm(u2)
tangent2 = u3/norm(u3)
push!(ntcoords, [normal tangent1 tangent2])
end
element["normal-tangential coordinates"] = ntcoords
end
src/lagrange.jl
+21 -1
View File
@@ -27,6 +27,10 @@ function calculate_lagrange_basis_coefficients(P, X)
return inv(A)'
end
function refcoords(X::Matrix)
return Vector{Float64}[X[:,i] for i=1:size(X,2)]
end
"""
Create new Lagrange element
@@ -37,7 +41,10 @@ Examples
macro create_lagrange_element(element_name, element_description, X, P)
eltype = esc(element_name)
quote
global get_basis, get_dbasis
global get_basis, get_dbasis,
get_reference_element_coordinates,
get_reference_element_midpoint
#basis, dbasis = calculate_lagrange_basis($P, $X)
C = calculate_lagrange_basis_coefficients($P, $X)
basis(xi) = C*$P(xi)
@@ -49,6 +56,15 @@ macro create_lagrange_element(element_name, element_description, X, P)
return basis(xi)'
end
XX = refcoords($X)
function get_reference_element_coordinates(::Type{$eltype})
return XX
end
XXX = vec(mean($X, 2))
function get_reference_element_midpoint(::Type{$eltype})
return XXX
end
#=
function get_dbasis(::Type{$eltype}, xi::Vector{Float64})
return dbasis(xi)'
@@ -81,6 +97,10 @@ end
0.0 0.0 1.0],
(xi) -> [1.0, xi[1], xi[2]])
#function get_reference_element_midpoint(::Type{Tri3})
# return [1.0/3.0, 1.0/3.0]
#end
@create_lagrange_element(Tri6, "6 node quadratic triangle element",
[0.0 1.0 0.0 0.5 0.5 0.0
0.0 0.0 1.0 0.0 0.5 0.5],
src/mortar.jl
+22 -8
View File
@@ -11,7 +11,7 @@ function project_from_slave_to_master{S,M}(slave::Element{S}, master::Element{M}
# slave side geometry and normal direction at xi1
X1 = slave("geometry", xi1, time)
N1 = slave("nodal ntsys", xi1, time)[:,1]
N1 = slave("normal-tangential coordinates", xi1, time)[:,1]
# master side geometry at xi2
#master_basis = master.basis.data.basis
@@ -59,7 +59,7 @@ function project_from_master_to_slave{S,M}(slave::Element{S}, master::Element{M}
# slave side geometry and normal direction at xi1
slave_geometry = slave("geometry")(time)
slave_normals = slave("nodal ntsys")(time)
slave_normals = slave("normal-tangential coordinates")(time)
#slave_basis = slave.basis.data.basis
#slave_dbasis = slave.basis.data.dbasis
slave_basis(xi) = get_basis(S, [xi])
@@ -86,7 +86,7 @@ function project_from_master_to_slave{S,M}(slave::Element{S}, master::Element{M}
end
#X1(xi1) = slave_basis("geometry", [xi1], time)
#N1(xi1) = slave_basis("nodal ntsys", [xi1], time)[:,1]
#N1(xi1) = slave_basis("normal-tangential coordinates", [xi1], time)[:,1]
#master_basis = get_basis(master)
@@ -187,11 +187,16 @@ Notes
"""
# function create_auxiliary_plane(x, ximp, normals, basis)
function create_auxiliary_plane(element::Element{Tri3}, time::Real)
proj(u, v) = dot(v, u) / dot(u, u) * u
xi = [1.0/3.0, 1.0/3.0]
function create_auxiliary_plane{E}(element::Element{E}, time::Real)
# proj(u, v) = dot(v, u) / dot(u, u) * u
# xi = [1.0/3.0, 1.0/3.0]
xi = get_reference_element_midpoint(E)
x0 = element("geometry", xi, time)
n = element("nodal ntsys", xi, time)[:, 1]
ntbasis = element("normal-tangential coordinates", xi, time)
return x0, ntbasis
#=
n = element("normal-tangential coordinates", xi, time)[:, 1]
n /= norm(n)
# gram-schmidt
u1 = n
@@ -204,6 +209,7 @@ function create_auxiliary_plane(element::Element{Tri3}, time::Real)
t2 = u3/norm(u3)
new_basis = [n t1 t2]
return x0, new_basis
=#
end
"""
@@ -243,7 +249,14 @@ function project_point_to_auxiliary_plane(p::Vector, x0::Vector, Q::Matrix)
n = Q[:,1]
ph = p - dot(p-x0, n)*n
qproj = Q'*(ph-x0)
@assert isapprox(qproj[1], 0.0)
if !isapprox(qproj[1], 0.0; atol=1.0e-12)
info("project_point_to_auxiliary_plane(): point not projected correctly.")
info("p: $p")
info("x0: $x0")
info("Q: \n$Q")
info("qproj: $qproj")
error("Failed to project point to auxiliary plane.")
end
return qproj[2:3]
end
@@ -626,6 +639,7 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly, problem::BoundaryPro
N2 = master_element(xi_projected, time)
S = w*N1'*N1
M = w*(N1'*N2)'
# M = w*N1'*N2
# FIXME: why this needs now to be transpose?
# assembly / repeat
for i=1:field_dim
test/test_elements.jl
+17 -1
View File
@@ -6,7 +6,8 @@ module ElementTests
using JuliaFEM.Test
using JuliaFEM.Core: AbstractElement, Element, Field, FieldSet, test_element
import JuliaFEM.Core: get_basis, get_dbasis
using JuliaFEM.Core: Tri3
import JuliaFEM.Core: get_basis, get_dbasis, calculate_normal_tangential_coordinates!
import Base: size
""" Prototype element
@@ -75,4 +76,19 @@ function test_interpolate()
# @test isapprox(gradT, 1/2*gradT_expected)
end
function test_calculate_normal_tangential_coordinates()
el = Tri3([1, 2, 3])
el["geometry"] = Vector{Float64}[
[0.0, 0.0, 0.0],
[1.0, 0.0, 0.0],
[0.0, 1.0, 0.0]]
calculate_normal_tangential_coordinates!(el, 0.0)
n = [0.0 0.0 1.0]'
t1 = [1.0 0.0 0.0]'
t2 = [0.0 1.0 0.0]'
R = [n t1 t2]
@test isapprox(el("normal-tangential coordinates", [0.0, 0.0], 0.0), R)
end
#test_calculate_normal_tangential_coordinates()
end
test/test_mortar.jl
+66 -25
View File
@@ -15,7 +15,8 @@ using JuliaFEM.Core: project_from_slave_to_master, project_from_master_to_slave
using JuliaFEM.Core: create_auxiliary_plane, project_point_to_auxiliary_plane,
get_edge_intersections, get_points_inside_triangle,
clip_polygon, calculate_polygon_centerpoint,
project_point_from_plane_to_surface, assemble
project_point_from_plane_to_surface, assemble,
calculate_normal_tangential_coordinates!
function get_test_2d_model()
@@ -42,13 +43,13 @@ function get_test_2d_model()
slave1 = Seg2([10, 11])
slave1["geometry"] = Vector[N[10], N[11]]
# should be n = [0 -1]' and t = [1 0]'
slave1["nodal ntsys"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
slave1["normal-tangential coordinates"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
slave1["master elements"] = Element[master1, master2]
slave2 = Seg2([11, 12])
slave2["geometry"] = Vector[N[11], N[12]]
# should be n = [0 -1]' and t = [1 0]'
slave2["nodal ntsys"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
slave2["normal-tangential coordinates"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
slave2["master elements"] = Element[master1, master2]
return [slave1, slave2], [master1, master2]
@@ -88,7 +89,7 @@ function test_calc_flat_2d_projection_rotated()
master1["geometry"] = Vector{Float64}[[0.0, 1.0], [0.0, 0.0]]
slave1 = Seg2([1, 2])
slave1["geometry"] = Vector{Float64}[[0.0, 0.0], [0.0, 1.0]]
slave1["nodal ntsys"] = Matrix{Float64}[[1.0 0.0; 0.0 1.0], [1.0 0.0; 0.0 1.0]]
slave1["normal-tangential coordinates"] = Matrix{Float64}[[1.0 0.0; 0.0 1.0], [1.0 0.0; 0.0 1.0]]
xi = project_from_master_to_slave(slave1, master1, [-1.0])
info("xi = $xi")
@test xi == [ 1.0]
@@ -209,7 +210,7 @@ function test_2d_mortar_multiple_bodies_multiple_dirichlet_bc()
slave1 = Seg2([5, 6])
slave1["geometry"] = Vector[N[5], N[6]]
slave1["nodal ntsys"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
slave1["normal-tangential coordinates"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
slave1["master elements"] = Element[master1]
boundary3 = MortarProblem("displacement", 2)
@@ -237,13 +238,19 @@ end
function test_2d_mortar_three_bodies_shared_nodes()
N = Vector[
[0.0, 0.0], [2.0, 0.0],
[0.0, 1.0], [2.0, 1.0],
[0.0, 1.0], [1.0, 1.0],
[0.0, 2.0], [1.0, 2.0],
[1.0, 1.0], [2.0, 1.0],
[1.0, 2.0], [2.0, 2.0]]
N = Dict{Int, Vector{Float64}}(
1 => [0.0, 0.0],
2 => [2.0, 0.0],
3 => [0.0, 1.0],
4 => [2.0, 1.0],
5 => [0.0, 1.0],
6 => [1.3, 1.0],
7 => [0.0, 2.0],
8 => [1.3, 2.0],
9 => [1.3, 1.0],
10 => [2.0, 1.0],
11 => [1.3, 2.0],
12 => [2.0, 2.0])
e1 = Quad4([1, 2, 4, 3])
e1["geometry"] = Vector[N[1], N[2], N[4], N[3]]
@@ -307,7 +314,7 @@ function test_2d_mortar_three_bodies_shared_nodes()
slave1 = Seg2([5, 6])
slave1["geometry"] = Vector[N[5], N[6]]
slave1["nodal ntsys"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
slave1["normal-tangential coordinates"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
slave1["master elements"] = Element[master1]
bc3 = MortarProblem("displacement", 2)
push!(bc3, slave1)
@@ -315,7 +322,7 @@ function test_2d_mortar_three_bodies_shared_nodes()
# mortar boundary between body 1 and body 3
slave2 = Seg2([9, 10])
slave2["geometry"] = Vector[N[9], N[10]]
slave2["nodal ntsys"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
slave2["normal-tangential coordinates"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
slave2["master elements"] = Element[master1]
bc4 = MortarProblem("displacement", 2)
push!(bc4, slave2)
@@ -326,8 +333,8 @@ function test_2d_mortar_three_bodies_shared_nodes()
slave3 = Seg2([6, 8])
slave3["geometry"] = Vector[N[6], N[8]]
#slave3["nodal ntsys"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
slave3["nodal ntsys"] = Matrix[rotation_matrix(0.0), rotation_matrix(0.0)]
#slave3["normal-tangential coordinates"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
slave3["normal-tangential coordinates"] = Matrix[rotation_matrix(0.0), rotation_matrix(0.0)]
slave3["master elements"] = Element[master2]
bc5 = MortarProblem("displacement", 2)
push!(bc5, slave3)
@@ -346,12 +353,15 @@ function test_2d_mortar_three_bodies_shared_nodes()
# launch solver
solver.method = :UMFPACK
solver.name = "test_2d_mortar_three_bodies_shared_nodes"
solver.dump_matrices = true
call(solver, 0.0)
disp = e2("displacement", [1.0, 1.0], 0.0)
info("displacement at tip: $disp")
X = e3("geometry", [1.0, 1.0], 0.0)
u = e3("displacement", [1.0, 1.0], 0.0)
info("displacement at $X: $u")
# code aster verification, two_elements.comm
@test isapprox(disp, [3.17431158889468E-02, -2.77183037855653E-01])
@test isapprox(u, [2*3.17431158889468E-02, -2.77183037855653E-01])
end
#test_2d_mortar_three_bodies_shared_nodes()
@@ -368,7 +378,7 @@ function test_auxiliary_plane_transforms()
0.0 0.0 1.0
1.0 0.0 0.0]
e1["geometry"] = Vector{Float64}[nodes[1], nodes[2], nodes[3]]
e1["nodal ntsys"] = Matrix{Float64}[R, R, R]
e1["normal-tangential coordinates"] = Matrix{Float64}[R, R, R]
time::Real = 0.0
x0, Q = create_auxiliary_plane(e1, time)
info("x0 = $x0")
@@ -466,6 +476,7 @@ end
#test_calculate_polygon_centerpoint()
function test_assemble_3d_problem()
nodes = Vector{Float64}[
[0.0, 0.0, 0.0],
@@ -478,12 +489,14 @@ function test_assemble_3d_problem()
mel["geometry"] = Vector{Float64}[nodes[4], nodes[5], nodes[6]]
sel = Tri3([1, 2, 3])
sel["geometry"] = Vector{Float64}[nodes[1], nodes[2], nodes[3]]
R = [0.0 1.0 0.0
0.0 0.0 1.0
1.0 0.0 0.0]
sel["nodal ntsys"] = Matrix{Float64}[R, R, R]
# Rv = [0.0 1.0 0.0
# 0.0 0.0 1.0
# 1.0 0.0 0.0]
# sel["normal-tangential coordinates"] = Matrix{Float64}[Rv, Rv, Rv]
calculate_normal_tangential_coordinates!(sel, 0.0)
sel["master elements"] = Element[mel]
prob = MortarProblem("temperature", 1)
push!(prob, sel)
stiffness_matrix = full(assemble(prob, 0.0).stiffness_matrix)
info("stiffness matrix for this problem:\n$stiffness_matrix")
@@ -491,8 +504,36 @@ function test_assemble_3d_problem()
B = [D -M] # slave dofs are first in this.
info("expected matrix for this problem:\n$B")
@test isapprox(stiffness_matrix, B)
# rotate and translate surface and check that we are still having same results
Rx(t) = [
1.0 0.0 0.0
0.0 cos(t) -sin(t)
0.0 sin(t) cos(t)]
Ry(t) = [
cos(t) 0.0 sin(t)
0.0 1.0 0.0
-sin(t) 0.0 cos(t)
]
Rz(t) = [
cos(t) -sin(t) 0.0
sin(t) cos(t) 0.0
0.0 0.0 1.0]
T = [1.0, 1.0, 1.0]
tx = pi/3.0
ty = pi/4.0
tz = pi/5.0
for node in nodes
node[:] = Rz(tz)*Ry(ty)*Rx(tx)*node + T
end
calculate_normal_tangential_coordinates!(sel, 0.0)
stiffness_matrix = full(assemble(prob, 0.0).stiffness_matrix)
info("sel midpnt: ", sel("geometry", [1/3, 1/3], 0.0))
info("nt basis: ", sel("normal-tangential coordinates", [1/3, 1/3], 0.0))
@test isapprox(stiffness_matrix, B)
end
#test_assemble_3d_problem()
test_assemble_3d_problem()
end