Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-09-25 19:36:58 +00:00
Code Issues Packages Projects Releases Wiki Activity

use FEMBase v0.1.x (#185)

Browse Source 
Lots of stuff moved from JuliaFEM.jl to FEMBase.jl.
This commit is contained in:
Jukka Aho
2018-01-19 21:01:33 +07:00
committed by GitHub
parent 35307b12ad
commit d7bef419ed
45 changed files with 232 additions and 1246 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/JuliaFEM.jl
+1 -14
View File
@@ -9,18 +9,6 @@ This is JuliaFEM -- Finite Element Package
module JuliaFEM
using FEMBase
using FEMBase: SparseMatrixCOO, SparseVectorCOO, Node, BasisInfo,
Discrete, Variable, TimeVariant, TimeInvariant, Field,
DCTI, DVTI, DCTV, DVTV, CCTI, CVTI, CCTV, CVTV, Increment,
IP, AbstractProblem, IntegrationPoint
using FEMBase: is_field_problem, is_boundary_problem, get_elements,
get_connectivity, assemble_prehook!, assemble_posthook!,
get_parent_field_name, get_reference_coordinates,
get_assembly, get_nonzero_rows, get_nonzero_columns,
eval_basis!, get_basis, get_dbasis, grad!, get_dualbasis,
assemble_mass_matrix!, get_local_coordinates, inside,
get_element_type, filter_by_element_type, get_element_id,
optimize!, resize_sparse, resize_sparsevec
import FEMBase: get_unknown_field_name, get_unknown_field_dimension,
assemble!, update!, initialize!
@@ -107,7 +95,7 @@ end
include("deprecations.jl")
export SparseMatrixCOO, SparseVectorCOO, optimize!, resize_sparse
export Field, DCTI, DVTI, DCTV, DVTV, CCTI, CVTI, CCTV, CVTV, Increment
export DCTI, DVTI, DCTV, DVTV, CCTI, CVTI, CCTV, CVTV, Increment
export FieldProblem, BoundaryProblem, Problem, Node, Element, Assembly
export Poi1, Seg2, Seg3, Tri3, Tri6, Tri7, Quad4, Quad8, Quad9,
Tet4, Tet10, Pyr5, Wedge6, Wedge15, Hex8, Hex20, Hex27
@@ -115,7 +103,6 @@ export update!, add_elements!, get_unknown_field_name, add!,
is_field_problem, is_boundary_problem, get_gdofs,
initialize!, get_integration_points, group_by_element_type,
get_unknown_field_dimension, get_connectivity
export get_nonzero_rows, get_local_coordinates, inside, IP, get_element_type,
get_elements, AbstractProblem, IntegrationPoint, filter_by_element_type,
get_element_id, get_nonzero_columns, resize_sparse, resize_sparsevec
src/io.jl
+1 -1
View File
@@ -429,7 +429,7 @@ function update_xdmf!(xdmf::Xdmf, problem::Problem, time::Float64, fields::Vecto
info("Xdmf: Saving topology of $nelements elements total, $nelement_types different element types.")
for element_type in element_types
elements = filter_by_element_type(element_type, all_elements)
elements = collect(filter_by_element_type(element_type, all_elements))
nelements = length(elements)
info("Xdmf: $nelements elements of type $element_type")
sort!(elements, by=get_element_id)
src/postprocess_utils.jl
+1 -1
View File
@@ -69,7 +69,7 @@ Return node ids + vector of values
function get_nodal_vector(elements::Vector, field_name::AbstractString, time::Float64)
f = Dict()
for element in elements
for (c, v) in zip(get_connectivity(element), element[field_name](time))
for (c, v) in zip(get_connectivity(element), element(field_name, time))
if haskey(f, c)
@assert isapprox(f[c], v)
end
src/problems_contact_2d.jl
+11 -11
View File
@@ -80,7 +80,7 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{1}}, ::T
la1 = slave_element("lambda", time)
n1 = slave_element("normal", time)
t1 = slave_element("tangent", time)
x1 = X1 + u1
x1 = map(+, X1, u1)
contact_area = 0.0
contact_error = 0.0
@@ -116,7 +116,7 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{1}}, ::T
nm = length(master_element)
X2 = master_element("geometry", time)
u2 = master_element("displacement", time)
x2 = X2 + u2
x2 = map(+, X2, u2)
# 3.3. loop integration points of one integration segment and calculate
# local mortar matrices
@@ -136,20 +136,20 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{1}}, ::T
Phi = Ae*N1
# project gauss point from slave element to master element in direction n_s
X_s = N1*X1 # coordinate in gauss point
n_s = N1*n1 # normal direction in gauss point
t_s = N1*t1 # tangent condition in gauss point
X_s = interpolate(N1, X1) # coordinate in gauss point
n_s = interpolate(N1, n1) # normal direction in gauss point
t_s = interpolate(N1, t1) # tangent condition in gauss point
n_s /= norm(n_s)
t_s /= norm(t_s)
xi_m = project_from_slave_to_master(master_element, X_s, n_s, time)
N2 = vec(get_basis(master_element, xi_m, time))
X_m = N2*X2
X_m = interpolate(N2, X2)
u_s = N1*u1
u_m = N2*u2
x_s = X_s + u_s
x_m = X_m + u_m
la_s = Phi*la1
u_s = interpolate(N1, u1)
u_m = interpolate(N2, u2)
x_s = map(+, X_s, u_s)
x_m = map(+, X_m, u_m)
la_s = interpolate(Phi, la1)
# virtual work
De += w*Phi*N1'
src/problems_contact_2d_autodiff.jl
+33 -28
View File
@@ -19,14 +19,20 @@ xi
projected master
"""
function project_from_master_to_slave{E<:MortarElements2D}(
function project_from_master_to_slave_ad{E<:MortarElements2D}(
slave_element::Element{E}, x1_::DVTI, n1_::DVTI, x2::Vector;
tol=1.0e-10, max_iterations=20, debug=false)
x1(xi1) = vec(get_basis(slave_element, [xi1], time))*x1_
dx1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*x1_
n1(xi1) = vec(get_basis(slave_element, [xi1], time))*n1_
dn1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*n1_
""" Multiply basis / dbasis at `xi` with field. """
function mul(func, xi, field)
B = func(slave_element, [xi], time)
return sum(B[i]*field[i] for i=1:length(B))
end
x1(xi1) = mul(get_basis, xi1, x1_)
dx1(xi1) = mul(get_dbasis, xi1, x1_)
n1(xi1) = mul(get_basis, xi1, n1_)
dn1(xi1) = mul(get_dbasis, xi1, n1_)
cross2(a, b) = cross([a; 0], [b; 0])[3]
R(xi1) = cross2(x1(xi1)-x2, n1(xi1))
dR(xi1) = cross2(dx1(xi1), n1(xi1)) + cross2(x1(xi1)-x2, dn1(xi1))
@@ -62,12 +68,12 @@ function project_from_master_to_slave{E<:MortarElements2D}(
end
function project_from_slave_to_master{E<:MortarElements2D}(
master_element::Element{E}, x1::Vector, n1::Vector, x2_::DVTI;
function project_from_slave_to_master_ad{E<:MortarElements2D}(
master_element::Element{E}, x1, n1, x2_;
tol=1.0e-10, max_iterations=20)
x2(xi2) = vec(get_basis(master_element, [xi2], time))*x2_
dx2(xi2) = vec(get_dbasis(master_element, [xi2], time))*x2_
x2(xi2) = interpolate(vec(get_basis(master_element, [xi2], time)), x2_)
dx2(xi2) = interpolate(vec(get_dbasis(master_element, [xi2], time)), x2_)
cross2(a, b) = cross([a; 0], [b; 0])[3]
R(xi2) = cross2(x2(xi2)-x1, n1)
dR(xi2) = cross2(dx2(xi2), n1)
@@ -119,8 +125,7 @@ function assemble!(problem::Problem{Contact}, time::Float64,
push!(S, conn...)
gdofs = get_gdofs(element, field_dim)
X_el = element("geometry", time)
u_el = Field(Vector[u[:,i] for i in conn])
x_el = X_el + u_el
x_el = tuple( (X_el[i] + u[:,j] for (i,j) in enumerate(conn))... )
#=
for ip in get_integration_points(element, 3)
dN = get_dbasis(element, ip, time)
@@ -157,10 +162,10 @@ function assemble!(problem::Problem{Contact}, time::Float64,
slave_element_nodes = get_connectivity(slave_element)
X1 = slave_element("geometry", time)
u1 = Field(Vector[u[:,i] for i in slave_element_nodes])
x1 = X1 + u1
la1 = Field(Vector[la[:,i] for i in slave_element_nodes])
n1 = Field(Vector[normals[:,i] for i in slave_element_nodes])
u1 = ((u[:,i] for i in slave_element_nodes)...)
x1 = ((Xi+ui for (Xi,ui) in zip(X1,u1))...)
la1 = ((la[:,i] for i in slave_element_nodes)...)
n1 = ((normals[:,i] for i in slave_element_nodes)...)
nnodes = size(slave_element, 2)
# construct dual basis
@@ -170,12 +175,12 @@ function assemble!(problem::Problem{Contact}, time::Float64,
master_element_nodes = get_connectivity(master_element)
X2 = master_element("geometry", time)
u2 = Field(Vector[u[:,i] for i in master_element_nodes])
x2 = X2 + u2
u2 = ((u[:,i] for i in master_element_nodes)...)
x2 = ((Xi+ui for (Xi,ui) in zip(X2,u2))...)
# calculate segmentation: we care only about endpoints
xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1])
xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[2])
xi1a = project_from_master_to_slave_ad(slave_element, field(x1), field(n1), x2[1])
xi1b = project_from_master_to_slave_ad(slave_element, field(x1), field(n1), x2[2])
xi1 = clamp.([xi1a; xi1b], -1.0, 1.0)
l = 1/2*abs(xi1[2]-xi1[1])
isapprox(l, 0.0) && continue # no contribution in this master element
@@ -200,8 +205,8 @@ function assemble!(problem::Problem{Contact}, time::Float64,
master_element_nodes = get_connectivity(master_element)
X2 = master_element("geometry", time)
u2 = Field(Vector[u[:,i] for i in master_element_nodes])
x2 = X2 + u2
u2 = ((u[:,i] for i in master_element_nodes)...)
x2 = ((Xi+ui for (Xi,ui) in zip(X2,u2))...)
#x1_midpoint = 1/2*(x1[1]+x1[2])
#x2_midpoint = 1/2*(x2[1]+x2[2])
@@ -209,8 +214,8 @@ function assemble!(problem::Problem{Contact}, time::Float64,
#distance > props.maximum_distance && continue
# calculate segmentation: we care only about endpoints
xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1])
xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[2])
xi1a = project_from_master_to_slave_ad(slave_element, field(x1), field(n1), x2[1])
xi1b = project_from_master_to_slave_ad(slave_element, field(x1), field(n1), x2[2])
xi1 = clamp.([xi1a; xi1b], -1.0, 1.0)
l = 1/2*abs(xi1[2]-xi1[1])
isapprox(l, 0.0) && continue # no contribution in this master element
@@ -229,15 +234,15 @@ function assemble!(problem::Problem{Contact}, time::Float64,
xi = ip.coords[1]
xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
N1 = vec(get_basis(slave_element, xi_s, time))
x_s = N1*x1 # coordinate in gauss point
n_s = N1*n1 # normal direction in gauss point
x_s = interpolate(N1, x1) # coordinate in gauss point
n_s = interpolate(N1, n1) # normal direction in gauss point
t_s = Q'*n_s # tangent direction in gauss point
xi_m = project_from_slave_to_master(master_element, x_s, n_s, x2)
xi_m = project_from_slave_to_master_ad(master_element, x_s, n_s, x2)
N2 = vec(get_basis(master_element, xi_m, time))
x_m = N2*x2
x_m = interpolate(N2, x2)
Phi = Ae*N1
la_s = Phi*la1 # traction force in gauss point
la_s = interpolate(Phi, la1) # traction force in gauss point
gn = -dot(n_s, x_s - x_m) # normal gap
fc[:,slave_element_nodes] += w*la_s*N1'
src/problems_contact_3d.jl
+18 -18
View File
@@ -83,13 +83,13 @@ function create_contact_segmentation(slave_element, master_elements, x0, n0, tim
result = []
x1 = slave_element("geometry", time)
if deformed
x1 += slave_element("displacement", time)
x1 = map(+, x1, slave_element("displacement", time))
end
S = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in x1]
for master_element in master_elements
x2 = master_element("geometry", time)
if deformed
x2 += master_element("displacement", time)
x2 = map(+, x2, master_element("displacement", time))
end
M = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in x2]
P = get_polygon_clip(S, M, n0)
@@ -115,7 +115,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time
nsl = length(slave_element)
X1 = slave_element("geometry", time)
u1 = slave_element("displacement", time)
x1 = X1 + u1
x1 = map(+, X1, u1)
n1 = slave_element("normal", time)
la = slave_element("lambda", time)
@@ -124,8 +124,8 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time
# project slave nodes to auxiliary plane (x0, Q)
xi = get_mean_xi(slave_element)
N = vec(get_basis(slave_element, xi, time))
x0 = N*X1
n0 = N*n1
x0 = interpolate(N, X1)
n0 = interpolate(N, n1)
# create contact segmentation
segmentation = create_contact_segmentation(slave_element, slave_element("master elements", time), x0, n0, time)
@@ -147,7 +147,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time
# loop integration cells
for cell in get_cells(P, C0)
virtual_element = Element(Tri3, Int[])
update!(virtual_element, "geometry", cell)
update!(virtual_element, "geometry", tuple(cell...))
for ip in get_integration_points(virtual_element, 3)
detJ = virtual_element(ip, time, Val{:detJ})
w = ip.weight*detJ
@@ -173,7 +173,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time
nm = length(master_element)
X2 = master_element("geometry", time)
u2 = master_element("displacement", time)
x2 = X2 + u2
x2 = map(+, X2, u2)
De = zeros(nsl, nsl)
Me = zeros(nsl, nm)
@@ -183,7 +183,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time
# loop integration cells
for cell in get_cells(P, C0)
virtual_element = Element(Tri3, Int[])
update!(virtual_element, "geometry", cell)
update!(virtual_element, "geometry", tuple(cell...))
# loop integration point of integration cell
for ip in get_integration_points(virtual_element, 3)
@@ -202,8 +202,8 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time
De += w*Phi*N1'
Me += w*Phi*N2'
x_s = N1*(X1+u1)
x_m = N2*(X2+u2)
x_s = interpolate(N1, map(+,X1,u1))
x_m = interpolate(N2, map(+,X2,u2))
ge += w*vec((x_m-x_s)*Phi')
end # integration points done
@@ -282,8 +282,8 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time
# create auxiliary plane
xi = get_mean_xi(sub_slave_element)
N = vec(get_basis(sub_slave_element, xi, time))
x0 = N*X1
n0 = N*n1
x0 = interpolate(N, X1)
n0 = interpolate(N, n1)
# project slave nodes to auxiliary plane
S = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X1]
@@ -323,7 +323,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time
# 4. loop integration cells
for cell in get_cells(P, C0)
virtual_element = Element(Tri3, Int[])
update!(virtual_element, "geometry", cell)
update!(virtual_element, "geometry", tuple(cell...))
for ip in get_integration_points(virtual_element, 3)
detJ = virtual_element(ip, time, Val{:detJ})
w = ip.weight*detJ
@@ -358,8 +358,8 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time
# create auxiliary plane
xi = get_mean_xi(sub_slave_element)
N = vec(get_basis(sub_slave_element, xi, time))
x0 = N*X1
n0 = N*n1
x0 = interpolate(N, X1)
n0 = interpolate(N, n1)
# project slave nodes to auxiliary plane
S = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X1]
@@ -406,7 +406,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time
# 4. loop integration cells
for cell in get_cells(P, C0)
virtual_element = Element(Tri3, Int[])
update!(virtual_element, "geometry", cell)
update!(virtual_element, "geometry", tuple(cell...))
# 5. loop integration point of integration cell
for ip in get_integration_points(virtual_element, 3)
@@ -429,8 +429,8 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time
us = slave_element("displacement", time)
um = master_element("displacement", time)
xs = N1*(Xs+us)
xm = N2*(Xs+um)
xs = interpolate(N1, map(+,Xs,us))
xm = interpolate(N2, map(+,Xs,um))
ge += w*vec((xm-xs)*Phi')
end # integration points done
src/problems_dirichlet.jl
+15
View File
@@ -15,6 +15,21 @@ function Dirichlet()
Dirichlet(:incremental, false, false, 1)
end
""" Return dual basis transformation matrix Ae. """
function get_dualbasis(element::Element, time::Float64, order=1)
nnodes = length(element)
De = zeros(nnodes, nnodes)
Me = zeros(nnodes, nnodes)
for ip in get_integration_points(element, order)
detJ = element(ip, time, Val{:detJ})
w = ip.weight*detJ
N = element(ip, time)
De += w*diagm(vec(N))
Me += w*N'*N
end
return De, Me, De*inv(Me)
end
function get_formulation_type(problem::Problem{Dirichlet})
return problem.properties.formulation
end
src/problems_heat.jl
+4 -4
View File
@@ -99,7 +99,7 @@ function assemble!{E<:Heat3DVolumeElements}(assembly::Assembly, problem::Problem
fq += w*N'*f
end
end
T = vec(element[field_name](time))
T = [interpolate(element[field_name], time)...]
fq -= K*T
add!(assembly.K, gdofs, gdofs, K)
add!(assembly.f, gdofs, fq)
@@ -136,7 +136,7 @@ function assemble!{E<:Heat3DSurfaceElements}(assembly::Assembly, problem::Proble
fq += w*N'*h*Tu
end
end
T = vec(element[field_name](time))
T = collect(element(field_name, time))
fq -= K*T
add!(assembly.K, gdofs, gdofs, K)
add!(assembly.f, gdofs, fq)
@@ -180,7 +180,7 @@ function assemble!{E<:Heat2DVolumeElements}(assembly::Assembly, problem::Problem
fq += w*N'*f
end
end
T = vec(element[field_name](time))
T = collect(element(field_name, time))
fq -= K*T
add!(assembly.K, gdofs, gdofs, K)
add!(assembly.f, gdofs, fq)
@@ -217,7 +217,7 @@ function assemble!{E<:Heat2DSurfaceElements}(assembly::Assembly, problem::Proble
fq += w*N'*h*Tu
end
end
T = vec(element[field_name](time))
T = collect(element(field_name, time))
fq -= K*T
add!(assembly.K, gdofs, gdofs, K)
add!(assembly.f, gdofs, fq)
src/problems_mortar.jl
+4 -4
View File
@@ -143,13 +143,13 @@ function diagnose_interface(problem::Problem{Mortar}, time::Float64)
info("Slave element connectivity = $slave_element_nodes")
nsl = length(slave_element)
X1 = slave_element("geometry", time)
n1 = Field([normals[j] for j in slave_element_nodes])
n1 = tuple(collect(normals[j] for j in slave_element_nodes)...)
# project slave nodes to auxiliary plane (x0, Q)
xi = get_mean_xi(slave_element)
N = vec(get_basis(slave_element, xi, time))
x0 = N*X1
n0 = N*n1
x0 = interpolate(N,X1)
n0 = interpolate(N,n1)
info("Auxiliary plane x0 = $x0, n0 = $n0")
S = Vector[project_vertex_to_auxiliary_plane(X1[i], x0, n0) for i=1:nsl]
check_orientation!(S, n0)
@@ -211,7 +211,7 @@ function diagnose_interface(problem::Problem{Mortar}, time::Float64)
for (cell_id, cell) in enumerate(all_cells)
C_area = 0.0
virtual_element = Element(Tri3, Int[])
update!(virtual_element, "geometry", cell)
update!(virtual_element, "geometry", tuple(cell...))
# 5. loop integration point of integration cell
for ip in get_integration_points(virtual_element, 3)
src/problems_mortar_2d.jl
+14 -14
View File
@@ -24,12 +24,12 @@ function get_slave_elements(problem::Problem)
end
function project_from_master_to_slave{E<:MortarElements2D}(slave_element::Element{E}, x2, time)
x1_ = slave_element["geometry"](time)
n1_ = slave_element["normal"](time)
x1(xi1) = vec(get_basis(slave_element, [xi1], time))*x1_
dx1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*x1_
n1(xi1) = vec(get_basis(slave_element, [xi1], time))*n1_
dn1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*n1_
x1_ = slave_element("geometry", time)
n1_ = slave_element("normal", time)
x1(xi1) = interpolate(vec(get_basis(slave_element, [xi1], time)), x1_)
dx1(xi1) = interpolate(vec(get_dbasis(slave_element, [xi1], time)), x1_)
n1(xi1) = interpolate(vec(get_basis(slave_element, [xi1], time)), n1_)
dn1(xi1) = interpolate(vec(get_dbasis(slave_element, [xi1], time)), n1_)
R(xi1) = cross2(x1(xi1)-x2, n1(xi1))
dR(xi1) = cross2(dx1(xi1), n1(xi1)) + cross2(x1(xi1)-x2, dn1(xi1))
xi1 = nothing
@@ -55,9 +55,9 @@ function project_from_master_to_slave{E<:MortarElements2D}(slave_element::Elemen
end
function project_from_slave_to_master{E<:MortarElements2D}(master_element::Element{E}, x1, n1, time)
x2_ = master_element["geometry"](time)
x2(xi2) = vec(get_basis(master_element, [xi2], time))*x2_
dx2(xi2) = vec(get_dbasis(master_element, [xi2], time))*x2_
x2_ = master_element("geometry", time)
x2(xi2) = interpolate(vec(get_basis(master_element, [xi2], time)), x2_)
dx2(xi2) = interpolate(vec(get_dbasis(master_element, [xi2], time)), x2_)
cross2(a, b) = cross([a; 0], [b; 0])[3]
R(xi2) = cross2(x2(xi2)-x1, n1)
dR(xi2) = cross2(dx2(xi2), n1)
@@ -173,11 +173,11 @@ function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Ty
N1 = vec(get_basis(slave_element, xi_s, time))
Phi = Ae*N1
# project gauss point from slave element to master element in direction n_s
X_s = N1*X1 # coordinate in gauss point
n_s = N1*n1 # normal direction in gauss point
X_s = interpolate(N1, X1) # coordinate in gauss point
n_s = interpolate(N1, n1) # normal direction in gauss point
xi_m = project_from_slave_to_master(master_element, X_s, n_s, time)
N2 = vec(get_basis(master_element, xi_m, time))
X_m = N2*X2
X_m = interpolate(N2, X2)
De += w*Phi*N1'
Me += w*Phi*N2'
if props.adjust
@@ -187,8 +187,8 @@ function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Ty
norm(mean(X1) - X2[2]) / norm(X1[2] - X1[1]) < props.distval || continue
u1 = slave_element("displacement", time)
u2 = master_element("displacement", time)
x_s = X_s + N1*u1
x_m = X_m + N2*u2
x_s = X_s + interpolate(N1, u1)
x_m = X_m + interpolate(N2, u2)
ge += w*vec((x_m-x_s)*Phi')
end
end
src/problems_mortar_2d_autodiff.jl
+30 -30
View File
@@ -5,14 +5,14 @@ using ForwardDiff
# forwarddiff version of mesh tying in 2d
function project_from_master_to_slave{E<:MortarElements2D}(
slave_element::Element{E}, x1_::DVTI, n1_::DVTI, x2::Vector, time::Float64;
function project_from_master_to_slave_ad{E<:MortarElements2D}(
slave_element::Element{E}, x1_, n1_, x2, time;
tol=1.0e-10, max_iterations=20)
x1(xi1) = vec(get_basis(slave_element, [xi1], time))*x1_
dx1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*x1_
n1(xi1) = vec(get_basis(slave_element, [xi1], time))*n1_
dn1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*n1_
x1(xi1) = interpolate(vec(get_basis(slave_element, [xi1], time)), x1_)
dx1(xi1) = interpolate(vec(get_dbasis(slave_element, [xi1], time)), x1_)
n1(xi1) = interpolate(vec(get_basis(slave_element, [xi1], time)), n1_)
dn1(xi1) = interpolate(vec(get_dbasis(slave_element, [xi1], time)), n1_)
cross2(a, b) = cross([a; 0], [b; 0])[3]
R(xi1) = cross2(x1(xi1)-x2, n1(xi1))
dR(xi1) = cross2(dx1(xi1), n1(xi1)) + cross2(x1(xi1)-x2, dn1(xi1))
@@ -37,12 +37,12 @@ function project_from_master_to_slave{E<:MortarElements2D}(
end
function project_from_slave_to_master{E<:MortarElements2D}(
master_element::Element{E}, x1::Vector, n1::Vector, x2_::DVTI, time::Float64;
function project_from_slave_to_master_ad{E<:MortarElements2D}(
master_element::Element{E}, x1, n1, x2_, time;
tol=1.0e-10, max_iterations=20)
x2(xi2) = vec(get_basis(master_element, [xi2], time))*x2_
dx2(xi2) = vec(get_dbasis(master_element, [xi2], time))*x2_
x2(xi2) = interpolate(vec(get_basis(master_element, [xi2], time)), x2_)
dx2(xi2) = interpolate(vec(get_dbasis(master_element, [xi2], time)), x2_)
cross2(a, b) = cross([a; 0], [b; 0])[3]
R(xi2) = cross2(x2(xi2)-x1, n1)
dR(xi2) = cross2(dx2(xi2), n1)
@@ -93,8 +93,8 @@ function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Ty
conn = get_connectivity(element)
push!(S, conn...)
X1 = element("geometry", time)
u1 = Field([u[:,i] for i in conn])
x1 = X1 + u1
u1 = ((u[:,i] for i in conn)...)
x1 = map(+, X1, u1)
dN = get_dbasis(element, [0.0], time)
tangent = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
for nid in conn
@@ -123,11 +123,11 @@ function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Ty
nsl = length(slave_element)
slave_element_nodes = get_connectivity(slave_element)
X1 = slave_element["geometry"](time)
u1 = Field(Vector[u[:,i] for i in slave_element_nodes])
x1 = X1 + u1
la1 = Field(Vector[la[:,i] for i in slave_element_nodes])
n1 = Field(Vector[normals[:,i] for i in slave_element_nodes])
X1 = slave_element("geometry", time)
u1 = ((u[:,i] for i in slave_element_nodes)...)
x1 = map(+, X1, u1)
la1 = ((la[:,i] for i in slave_element_nodes)...)
n1 = ((normals[:,i] for i in slave_element_nodes)...)
# 3. loop all master elements
@@ -136,12 +136,12 @@ function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Ty
nm = length(master_element)
master_element_nodes = get_connectivity(master_element)
X2 = master_element("geometry", time)
u2 = Field(Vector[u[:,i] for i in master_element_nodes])
x2 = X2 + u2
u2 = ((u[:,i] for i in master_element_nodes)...)
x2 = map(+, X2, u2)
# 3.1 calculate segmentation
xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1], time)
xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[2], time)
xi1a = project_from_master_to_slave_ad(slave_element, x1, n1, x2[1], time)
xi1b = project_from_master_to_slave_ad(slave_element, x1, n1, x2[2], time)
# xi1a = project_from_master_to_slave(slave_element, X2[1], time)
# xi1b = project_from_master_to_slave(slave_element, X2[2], time)
xi1 = clamp.([xi1a; xi1b], -1.0, 1.0)
@@ -181,20 +181,20 @@ function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Ty
N1 = vec(get_basis(slave_element, xi_s, time))
Phi = Ae*N1
# project gauss point from slave element to master element in direction n_s
x_s = N1*x1 # coordinate in gauss point
n_s = N1*n1 # normal direction in gauss point
x_s = interpolate(N1, x1) # coordinate in gauss point
n_s = interpolate(N1, n1) # normal direction in gauss point
#xi_m = project_from_slave_to_master(master_element, X_s, n_s, time)
xi_m = project_from_slave_to_master(master_element, x_s, n_s, x2, time)
xi_m = project_from_slave_to_master_ad(master_element, x_s, n_s, x2, time)
N2 = vec(get_basis(master_element, xi_m, time))
x_m = N2*x2
x_m = interpolate(N2, x2)
la_s = Phi*la1
la_s = interpolate(Phi, la1)
gn = dot(n_s, x_s-x_m)
u_s = N1*u1
u_m = N2*u2
X_s = N1*X1
X_m = N2*X2
u_s = interpolate(N1, u1)
u_m = interpolate(N2, u2)
X_s = interpolate(N1, X1)
X_m = interpolate(N2, X2)
fc[:,slave_element_nodes] += w*la_s*N1'
fc[:,master_element_nodes] -= w*la_s*N2'
src/problems_mortar_3d.jl
+29 -20
View File
@@ -154,14 +154,23 @@ function get_polygon_clip{T}(xs::Vector{T}, xm::Vector{T}, n::T)
end
""" Project some vertex p to surface of element E using Newton's iterations. """
function project_vertex_to_surface{E}(p::Vector, x0::Vector, n0::Vector,
element::Element{E}, x::DVTI, time::Real;
max_iterations::Int=10, iter_tol::Float64=1.0e-6)
function project_vertex_to_surface{E}(p, x0, n0,
element::Element{E}, x, time;
max_iterations=10, iter_tol=1.0e-6)
basis(xi) = get_basis(element, xi, time)
dbasis(xi) = get_dbasis(element, xi, time)
function dbasis(xi)
return get_dbasis(element, xi, time)
end
nnodes = length(element)
f(theta) = basis(theta[1:2])*x - theta[3]*n0 - p
L(theta) = inv3([dbasis(theta[1:2])*x -n0])
mul(a,b) = sum((a[:,i]*b[i]')' for i=1:length(b))
function f(theta)
b = [basis(theta[1:2])*collect(x)...;]
b = b - theta[3]*n0 - p
return b
end
L(theta) = inv3([mul(dbasis(theta[1:2]), x) -n0])
theta = zeros(3)
dtheta = zeros(3)
for i=1:max_iterations
@@ -371,8 +380,8 @@ function assemble!{E<:Union{Tri3, Quad4}}(problem::Problem{Mortar}, slave_elemen
xi = get_mean_xi(slave_element)
first_slave_element && debug("midpoint xi = $xi")
N = vec(get_basis(slave_element, xi, time))
x0 = N*X1
n0 = N*n1
x0 = interpolate(N, X1)
n0 = interpolate(N, n1)
S = Vector[project_vertex_to_auxiliary_plane(X1[i], x0, n0) for i=1:nsl]
master_elements = slave_element("master elements", time)
@@ -412,7 +421,7 @@ function assemble!{E<:Union{Tri3, Quad4}}(problem::Problem{Mortar}, slave_elemen
all_cells = get_cells(P, C0)
for cell in all_cells
virtual_element = Element(Tri3, Int[])
update!(virtual_element, "geometry", cell)
update!(virtual_element, "geometry", tuple(cell...))
for ip in get_integration_points(virtual_element, 3)
detJ = virtual_element(ip, time, Val{:detJ})
w = ip.weight*detJ
@@ -476,7 +485,7 @@ function assemble!{E<:Union{Tri3, Quad4}}(problem::Problem{Mortar}, slave_elemen
all_cells = get_cells(P, C0)
for cell in all_cells
virtual_element = Element(Tri3, Int[])
virtual_element.fields["geometry"] = DVTI(cell)
update!(virtual_element, "geometry", tuple(cell...))
# 5. loop integration point of integration cell
for ip in get_integration_points(virtual_element, 3)
@@ -499,8 +508,8 @@ function assemble!{E<:Union{Tri3, Quad4}}(problem::Problem{Mortar}, slave_elemen
if props.adjust && haskey(slave_element, "displacement") && haskey(master_element, "displacement")
u1 = slave_element("displacement", time)
u2 = master_element("displacement", time)
x_s = N1*(X1+u1)
x_m = N2*(X2+u2)
x_s = interpolate(N1, map(+,X1,u1))
x_m = interpolate(N2, map(+,X2,u2))
ge += w*vec((x_m-x_s)*Phi')
end
area += w
@@ -593,8 +602,8 @@ function assemble!{E<:Union{Tri6}}(problem::Problem{Mortar}, slave_element::Elem
xi = get_mean_xi(sub_slave_element)
first_slave_element && debug("midpoint xi = $xi")
N = vec(get_basis(sub_slave_element, xi, time))
x0 = N*X1
n0 = N*n1
x0 = interpolate(N, X1)
n0 = interpolate(N, n1)
# project slave nodes to auxiliary plane
S = Vector[project_vertex_to_auxiliary_plane(X1[i], x0, n0) for i=1:nsl]
@@ -633,7 +642,7 @@ function assemble!{E<:Union{Tri6}}(problem::Problem{Mortar}, slave_element::Elem
all_cells = get_cells(P, C0)
for cell in all_cells
virtual_element = Element(Tri3, Int[])
update!(virtual_element, "geometry", cell)
update!(virtual_element, "geometry", tuple(cell...))
for ip in get_integration_points(virtual_element, 3)
x_gauss = virtual_element("geometry", ip, time)
xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, Xs, time)
@@ -676,8 +685,8 @@ function assemble!{E<:Union{Tri6}}(problem::Problem{Mortar}, slave_element::Elem
xi = get_mean_xi(sub_slave_element)
first_slave_element && debug("midpoint xi = $xi")
N = vec(get_basis(sub_slave_element, xi, time))
x0 = N*X1
n0 = N*n1
x0 = interpolate(N, X1)
n0 = interpolate(N, n1)
# project slave nodes to auxiliary plane
S = Vector[project_vertex_to_auxiliary_plane(X1[i], x0, n0) for i=1:nsl]
@@ -736,7 +745,7 @@ function assemble!{E<:Union{Tri6}}(problem::Problem{Mortar}, slave_element::Elem
all_cells = get_cells(P, C0)
for cell in all_cells
virtual_element = Element(Tri3, Int[])
update!(virtual_element, "geometry", cell)
update!(virtual_element, "geometry", tuple(cell...))
# 5. loop integration point of integration cell
for ip in get_integration_points(virtual_element, 3)
@@ -758,8 +767,8 @@ function assemble!{E<:Union{Tri6}}(problem::Problem{Mortar}, slave_element::Elem
if props.adjust && haskey(slave_element, "displacement") && haskey(master_element, "displacement")
u1 = slave_element("displacement", time)
u2 = master_element("displacement", time)
xs = N1*(Xs+u1)
xm = N2*(Xm+u2)
xs = interpolate(N1, map(+,Xs,u1))
xm = interpolate(N2, map(+,Xm,u2))
ge += w*vec((xm-xs)*Phi')
end
area += w
src/solvers.jl
+4 -4
View File
@@ -14,7 +14,7 @@ type Solver{S<:AbstractSolver}
u :: Vector{Float64}
la :: Vector{Float64}
alpha :: Float64 # generalized alpha time integration coefficient
fields :: Dict{AbstractString, Field}
fields :: Dict{String, AbstractField}
properties :: S
end
@@ -314,9 +314,9 @@ function solve!(solver::Solver; empty_assemblies_before_solution=true, symmetric
end
if !haskey(solver, "fint")
solver.fields["fint"] = Field(time => f)
solver.fields["fint"] = field(solver.time => f)
else
update!(solver.fields["fint"], time => f)
update!(solver.fields["fint"], solver.time => f)
end
fint = solver.fields["fint"]
@@ -327,7 +327,7 @@ function solve!(solver::Solver; empty_assemblies_before_solution=true, symmetric
debug("Using generalized-α time integration, α=$alpha")
K = (1-alpha)*K
C1 = (1-alpha)*C1
f = (1-alpha)*f + alpha*fint[end-1].data
f = (1-alpha)*f + alpha*fint.data[end-1].second
end
ndofs = solver.ndofs
src/solvers_modal.jl
+2 -2
View File
@@ -327,7 +327,7 @@ function update_xdmf!(solver::Solver{Modal})
elcon_arrays = Dict()
@timeit "create topology arrays" for element_type in element_types
elements = filter_by_element_type(element_type, all_elements)
elements = collect(filter_by_element_type(element_type, all_elements))
nelements = length(elements)
eldim = length(element_type)
element_conn = zeros(Int, eldim, nelements)
@@ -384,7 +384,7 @@ function update_xdmf!(solver::Solver{Modal})
for element_type in element_types
timeit("save topology of element type $element_type") do
elements = filter_by_element_type(element_type, all_elements)
elements = collect(filter_by_element_type(element_type, all_elements))
nelements = length(elements)
element_ids = map(get_element_id, elements)
element_conn = elcon_arrays[element_type]