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

some mortar code

Browse Source
This commit is contained in:
Jukka Aho
2015-11-12 07:16:25 +02:00
parent 260a110269
commit 50a0018f63
6 changed files with 221 additions and 250 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CONTRIBUTING.rst
+30 -4
View File
@@ -78,15 +78,40 @@ Testing is made easy by using our `Makefile`. From there one founds convenient
functions `make test`, `make test_file` and `make test_function` to make testing
more rapid.
Git issues
----------
One possible workflow: See
`this post <http://superuser.com/questions/181517/how-to-execute-a-command-whenever-a-file-changes>_`:
.. code:: bash
while sleep_until_modified.sh src/equations.jl; do clear; make test_file FILE=test/test_equations.jl; done
or something similar. Every time file `src/equations.jl` is changed, tests from
file `test/test_equations.jl` are run.
Git spesific things
-------------------
Have done local changes, want to get latest updates and get "Cannot pull with
rebase: You have unstaged changes. Please commit or stash them."
http://stackoverflow.com/questions/23517464/error-cannot-pull-with-rebase-you-have-unstaged-changes
Solution: git stash + git pull + git pop
.. code-block:: bash
git stash
git pull
git stash pop
Local modifictions done, want to get clean file from repo
.. code-block:: bash
git checkout -- <filename>
Having untracked files not wanting to commit at all? Put them to `.gitignore` if
they are somehow generally unnecessary things. Or if they are "personal" (have
made e.g. notebooks not wanting to commit), add them to your private ignore rules
`.git/info/exclude` so that they don't show as untracked for you.
Use of UTF-8 characters in program code
---------------------------------------
@@ -140,5 +165,6 @@ Documentation
-------------
We use restructured text to document this project. Information how to write rst
format is described `here <http://sphinx-doc.org/rest.html>`_. See issue
`#49 <https://github.com/JuliaFEM/JuliaFEM.jl/issues/49>`_.
`#49 <https://github.com/JuliaFEM/JuliaFEM.jl/issues/49>`_. Keep line width max
80 characters.
docs/build_notebooks.jl
+13 -3
View File
@@ -63,24 +63,34 @@ function run_notebooks()
end
runtime = toc()
bn = "tutorials/$(ipynb[1:end-6])"
#try
# run(`ipython nbconvert tutorials/$ipynb --to rst --output=$bn`)
#catch error
# warn("unable to convert notebook to rst format")
# Base.showerror(Base.STDOUT, error)
#end
try
run(`ipython nbconvert tutorials/$ipynb --to rst --output=$bn`)
run(`ipython nbconvert tutorials/$ipynb --to html --output=$bn`)
catch error
warn("unable to convert notebook to rst format")
warn("unable to convert notebook to html format")
Base.showerror(Base.STDOUT, error)
end
try
run(`ipython nbconvert tutorials/$ipynb --to latex --output=$bn`)
catch error
warn("unable to convert notebook to tex format")
Base.showerror(Base.STDOUT, error)
end
try
run(`lualatex $bn.tex`)
run(`lualatex --output-directory=tutorials $bn.tex`)
catch error
warn("unable to convert notebook from tex to pdf")
Base.showerror(Base.STDOUT, error)
end
data = Dict("author" => "unknown", "status" => status, "runtime" => runtime,
"filename" => ipynb, "last_run" => time(), "description"=>"")
res = parse_rst("$bn.rst")
src/dirichlet.jl
+23 -10
View File
@@ -10,12 +10,29 @@ abstract DirichletEquation <: Equation
type DirichletProblem <: BoundaryProblem
unknown_field_name :: ASCIIString
unknown_field_dimension :: Int
equations :: Array{DirichletEquation, 1}
element_mapping :: Dict{DataType, DataType}
equations :: Vector{DirichletEquation}
element_mapping :: Dict{Element, Equation}
field_value :: Function
end
function DirichletProblem(dimension::Int, field_value::Function=(X)->[0.0,0.0,0.0])
""" Initialize new Dirichlet boundary condition.
Parameters
----------
dimension
dimension of unknown field
field_value
boundary function
Examples
--------
Create u(X) = 0.0 boundary condition for three-dimensional elasticity problem:
>>> u(X) = [0.0, 0.0, 0.0]
>>> bc = DirichletProblem(3, u)
"""
function DirichletProblem(dimension::Int=1, field_value::Function=(X)->[0.0,0.0,0.0])
element_mapping = nothing
if dimension == 1
element_mapping = Dict(
@@ -28,12 +45,10 @@ end
""" Dirichlet boundary condition element for 2 node line segment """
type DBC2D2 <: DirichletEquation
element :: Seg2
integration_points :: Array{IntegrationPoint, 1}
integration_points :: Vector{IntegrationPoint}
end
function DBC2D2(element::Seg2)
integration_points = [
IntegrationPoint([-sqrt(1/3)], 1.0),
IntegrationPoint([+sqrt(1/3)], 1.0)]
integration_points = default_integration_points(element)
if !haskey(element, "reaction force")
element["reaction force"] = zeros(1, 2)
end
@@ -41,9 +56,7 @@ function DBC2D2(element::Seg2)
end
Base.size(equation::DBC2D2) = (1, 2)
function calculate_local_assembly!(assembly::LocalAssembly, equation::DirichletEquation,
unknown_field_name::ASCIIString, time::Number=Inf,
problem=nothing)
function calculate_local_assembly!(assembly::LocalAssembly, equation::DirichletEquation, unknown_field_name::ASCIIString, time::Number=0.0, problem=nothing)
initialize_local_assembly!(assembly, equation)
element = get_element(equation)
basis = get_basis(element)
src/field_utils.jl
-140
View File
@@ -1,140 +0,0 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
# There are here for now until I figure a better place for them.
"""
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
src/mortar.jl
+72 -93
View File
@@ -1,101 +1,80 @@
"""
calculate "local" normals in elements, in a way that
n = Nᵢnᵢ gives some reasonable results for ξ ∈ [-1, 1]
"""
function calculate_normals!(el::Element, t, field_name=symbol("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
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
# Mortar projection integration
abstract MortarEquation <: Equation
"""
Alter normal field such that normals of adjacent elements are averaged.
Parameters
----------
node_csys
coordinate system in node, normal + tangent + "binormal"
element_pairs
m x s matrix of boolean values, indicating elements sharing
common surface. s is number of slave elements and m is number
of master elements.
"""
function average_normals!(elements, normal_field=symbol("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
type MortarProblem <: BoundaryProblem
unknown_field_name :: ASCIIString
unknown_field_dimension :: Int
equations :: Vector{MortarEquation}
element_mapping :: Dict{Element, MortarEquation}
master_elements :: Vector{Element} # mortar surface
node_csys :: Dict{Int, Matrix{Float64}}
element_pairs :: Matrix{Bool}
end
function MortarProblem(dimension::Int=1, equations=[], master_elements=[])
element_mapping = Dict(
Seg2 => MBC2D2,
)
MortarProblem("reaction force", dimension, equations, element_mapping, master_elements, Dict(), zeros(0,0))
end
""" Mortar boundary condition element for 2-dimensional problem, 2 node line segment. """
type MBC2D2 <: MortarEquation
element :: Seg2 # == non-mortar surface element
integration_points :: Vector{IntegrationPoint}
end
function MBC2D2(element::Seg2)
integration_points = default_integration_points(element)
if !haskey(element, "reaction force")
element["reaction force"] = zeros(1, 2)
end
MBC2D2(element, integration_points)
end
Base.size(equation::MBC2D2) = (1, 2)
function find_master_elements(slave_element, problem)
# find slave element "position" in element pairs matrix
all_elements = map((equation) -> get_element(equation), problem.equations)
seid = findfirst(slave_element, all_elements)
info("slave element id = $seid")
# find master element "positions" in element pairs matrix
meids = find(problem.element_pairs[:, seid])
info("master element ids = $meids")
# master elements
master_elements = problem.master_elements[meids]
return master_elements
end
function calculate_local_assembly!(assembly::LocalAssembly, equation::MortarEquation, unknown_field_name::ASCIIString, time::Number=0.0, problem=nothing)
# slave element = non-mortar element where integration happens
# master element = mortar element projected to non-mortar side
isa(problem, Void) && error("Cannot create projection without problem")
initialize_local_assembly!(assembly, equation)
slave_element = get_element(equation)
basis = get_basis(slave_element)
detJ = det(basis)
master_elements = find_master_elements(equation, problem)
for master_element in master_elements
for ip in get_integration_points(slave_element)
mortar_basis = 0 # ...
assembly.stiffness_matrix += w*basis'*basis
assembly.force_vector += w*N'*gn
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
""" Find projection from slave nodes to master element. """
function calc_projection_slave_nodes_to_master_element(sel, mel)
X1 = get_field(sel, :Geometry)
N1 = get_field(sel, :Normals)
X2(xi) = interpolate(mel, :Geometry, xi)
dX2(xi) = dinterpolate(mel, :Geometry, xi)
R(xi, k) = det([X2(xi) - X1[k] N1[k]]')
dR(xi, k) = det([dX2(xi) N1[k]]')
xi2 = Vector[[0.0], [0.0]]
for k=1:2
xi = xi2[k]
for i=1:3
dxi = -R(xi, k)/dR(xi, k)
xi += dxi
if abs(dxi) < 1.0e-9
break
end
end
xi2[k] = xi
end
clamp!(xi2, -1, 1)
return xi2
end
""" Find projection from master nodes to slave element. """
function calc_projection_master_nodes_to_slave_element(sel, mel)
X1(xi) = interpolate(sel, :Geometry, xi)
dX1(xi) = dinterpolate(sel, :Geometry, xi)
N1(xi) = interpolate(sel, :Normals, xi)
dN1(xi) = dinterpolate(sel, :Normals, xi)
X2 = get_field(mel, :Geometry)
R(xi, k) = det([X1(xi) - X2[k] N1(xi)]')
dR(xi, k) = det([dX1(xi) N1(xi)]') + det([X1(xi) - X2[k] dN1(xi)]')
xi1 = Vector[[0.0], [0.0]]
for k=1:2
xi = xi1[k]
for i=1:3
dxi = -R(xi, k)/dR(xi, k)
xi += dxi
if abs(dxi) < 1.0e-9
break
end
end
xi1[k] = xi
end
clamp!(xi1, -1, 1)
return xi1
end
function has_projection(sel, mel)
xi1 = calc_projection_master_nodes_to_slave_element(sel, mel)
l = abs(xi1[2]-xi1[1])[1]
return l > 1.0e-9
end
"""
Calculate projection between 1d boundary elements
"""
function calc_projection(sel, mel)
xi1 = calc_projection_master_nodes_to_slave_element(sel, mel)
xi2 = calc_projection_slave_nodes_to_master_element(sel, mel)
return xi1, xi2
end
test/test_mortar.jl
+83
View File
@@ -0,0 +1,83 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
module MortarTests
using JuliaFEM
using JuliaFEM.Test
function test_calc_flat_2d_assembly()
# this is hand calculated and given example in my thesis
N = Vector[
[0.0, 2.0], [1.0, 2.0], [2.0, 2.0],
[0.0, 0.0], [1.0, 0.0], [2.0, 0.0],
[0.0, 1.0], [5/4, 1.0], [2.0, 1.0],
[0.0, 1.0], [3/4, 1.0], [2.0, 1.0]]
slave1 = Seg2([10, 11])
slave1["geometry"] = Vector[N10, N11]
slave2 = Seg2([11, 12])
slave2["geometry"] = Vector[N11, N12]
master1 = Seg2([7, 8])
master1["geometry"] = Vector[N7, N8]
master2 = Seg2([8, 9])
master2["geometry"] = Vector[N8, N9]
problem = MortarProblem()
push!(problem, slave1)
push!(problem, slave2)
push!(problem.master_elements, master1)
push!(problem.master_elements, master2)
rotation_matrix(phi) = [cos(phi) -sin(phi); sin(phi) cos(phi)]
# should be n = [0 -1]' and t = [1 0]'
@test isapprox(rotation_matrix(-phi/2), [[0 -1]' [1 0]'])
problem.node_csys = Dict(
10 => rotation_matrix(-phi/2),
11 => rotation_matrix(-phi/2),
12 => rotation_matrix(-phi/2))
# first index = master element id
# second index = slave element id
problem.element_pairs = zeros(2, 2)
# first slave element connects to master element 1
problem.element_pairs[1, 1] = true
# second slave element connects to master element 1
problem.element_pairs[1, 2] = true
# second slave element connects to master element 2
problem.element_pairs[2, 2] = true
B_expected = zeros(12, 9)
S1 = [10, 11]
M1 = [7, 8]
B_expected[S1,S1] += [1/4 1/8; 1/8 1/4]
B_expected[S1,M1] += [3/10 3/40; 9/40 3/20]
la = initialize_local_assembly(problem)
calculate_local_assembly!(la, problem.equations[1], "reaction force", 0.0, problem=problem)
B = full(la.lhs)
@test isapprox(B, B_expected)
fill!(B_expected, 0.0)
S2 = [11, 12]
M2 = [7, 8]
B_expected[S2,S2] += [49/150 11/150; 11/150 2/75]
B_expected[S2,M2] += [13/150 47/150; 1/75 13/150]
S3 = [11, 12]
M3 = [8, 9]
B_expected[S3,S3] += [9/100 27/200; 27/200 39/100]
B_expected[S3,M3] += [3/20 3/40; 9/40 3/10]
la = initialize_local_assembly(problem)
calculate_local_assembly!(la, problem.equations[1], "reaction force", 0.0, problem=problem)
B = full(la.lhs)
@test isapprox(B, B_expected)
end
end