Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-09-22 18:52:16 +00:00
Code Issues Packages Projects Releases Wiki Activity

Bug/mortar discretization (#88)

Browse Source 
* new mortar segmentation tests which are failing

* test_problems_mortar_3d.jl: first test (Tet4) pass

* solvers.jl: diagonal of A is now properly filled, if that option is used. Another option is to remove zero rows from matrix system, which is on by default

* problems_mortar.jl: added new function diagnose_interface to calculate quantities from interface hopefully revealing bugs in calculation

* problems_mortar_3d.jl: added docstring for check_orientation! and removed flooding debug messages not helping to debug anything

* solvers.jl: Another way to solve Ax = b

* Refactored code to make implementation of Tri6 assemble! easier

* Patch test with linear Tet4 elements and quadratic Tet10 elements pass

When using quadratic elements, in polygon clipping algorithm element is divided to linear sub-elements as proposed in [Puso2008]. Interpolation of Lagrange multiplier space is done using quadratic shape functions.

References
----------

[Puso2008] Puso, Michael A., T. A. Laursen, and Jerome Solberg. "A segment-to-segment mortar contact method for quadratic elements and large deformations." Computer Methods in Applied Mechanics and Engineering 197.6 (2008): 555-566.

* increased coverage by adding diagnose_interface

* test using dual basis, failing for unknown reason

* Fixed dual basis construction for Mortar/Tet4

The coefficient matrix Ae for one particular slave element e is the result performing numerical integration on *all* integration cells associated with this element [Popp2013]. Ae cannot be calculated "cell-wise" like it was done before. Now patch test will pass also using `interface.properties.dual_basis = true` option. Partially integrated slave elements are supported as well.

References
----------

[Popp2013] Popp, Alexander, et al. "Improved robustness and consistency of 3D contact algorithms based on a dual mortar approach." Computer Methods in Applied Mechanics and Engineering 264 (2013): 67-80.

* Minor modifications to preprocess.jl

- removed two functions which are unimplemented (but maybe planned in future)
- added function create_node_set_from_element_set!, which can be used, like name suggests, to create a node set from nodes belonging to some set of elements.

* solvers.jl: now prints a list of overconstrained nodes which can be easily copy-pasted to problem.assembly.removed_dofs list to solver overconstrained situation manually

* Increase code coverage

Added a new test which tests dual basis 3d mortar + adjust option when using Tet4 in elasticity problem.

* Tet10 + Dual basis still failing, others are working

* mortar 3d low level tests

* linear surface element projection tests pass

* Introduced basis transform constant alpha

Tet10 + dual basis patch test still failing, but single element low level routine tests gives expected results with alpha=0.2

* added new integration rule FPG12 for triangular elements

* added drop_tolerance option to remove very small values from constraint matrices

* Introduced a basis transform matrix T

Constructing bi-orthogonal basis for quadratic surfaces is ill-conditioned. By doing a basis transform N' = N*T for slave side displacement vector it's possible to construct a bi-orthogonal basis in a same way than with linear elements. Setting alpha=0.2 ensures that quadratic basis functions are strictly positive in practical cases.

* fix 3d clipping test routine, accepts only 3d vertices

* dropped number of integration poitns from 12 to 7 in quadratic mortar surfaces intrestingly gives more accurate results, maybe something numerical error in FPG12 integration rule..?

* added two displacement patch tests + output writing for all cases

* %s/Int64/Int/g

* Changed test data location

* Fine tuning of logging levels
This commit is contained in:
Jukka Aho
2017-02-25 18:40:14 +02:00
committed by Tero Frondelius
parent a5c093c1d6
commit f275ce3767
10 changed files with 2026 additions and 264 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/integrate.jl
+48 -15
View File
@@ -148,23 +148,56 @@ function get_integration_points(element::TriangularElement, ::Type{Val{4}})
return zip(weights, points)
end
""" 7 point integration rule for triangular elements.
References
----------
Code Aster documentation, http://code-aster.org/doc/default/fr/man_r/r3/r3.01.01.pdf
"""
function get_integration_points(element::TriangularElement, ::Type{Val{5}})
weights = 0.5*[
0.22500000000000,
0.13239415278851,
0.13239415278851,
0.13239415278851,
0.12593918054483,
0.12593918054483,
0.12593918054483]
A = 0.470142064105115
B = 0.101286507323456
P1 = 0.066197076394253
P2 = 0.062969590272413
weights = [9/80, P1, P1, P1, P2, P2, P2]
points = Vector{Float64}[
[0.33333333333333, 0.33333333333333],
[0.47014206410511, 0.47014206410511],
[0.47014206410511, 0.05971587178977],
[0.05971587178977, 0.47014206410511],
[0.10128650732346, 0.10128650732346],
[0.10128650732346, 0.79742698535309],
[0.79742698535309, 0.10128650732346]]
[1/3, 1/3],
[A, A],
[1-2A, A],
[A, 1-2A],
[B, B],
[1-2B, B],
[B, 1-2B]]
return zip(weights, points)
end
""" 12 point integration fule for triangular elements.
References
----------
Code Aster documentation, http://code-aster.org/doc/default/fr/man_r/r3/r3.01.01.pdf
"""
function get_integration_points{E<:TriangularElement}(element::Element{E}, ::Type{Val{:FPG12}})
A = 0.063089014491502
B = 0.249286745170910
C = 0.310352451033785
D = 0.053145049844816
P1 = 0.025422453185103
P2 = 0.058393137863189
P3 = 0.041425537809187
weights = [P1, P1, P1, P2, P2, P2, P3, P3, P3, P3, P3, P3]
points = Vector{Float64}[
[A, A],
[1-2A, A],
[A, 1-2A],
[B, B],
[1-2B, B],
[B, 1-2B],
[C, D],
[D, C],
[1-C-D, C],
[1-C,D, D],
[C, 1-C-D],
[D, 1-C-D]]
return zip(weights, points)
end
src/preprocess.jl
+27 -31
View File
@@ -17,13 +17,13 @@ import Base: copy
using JuliaFEM
type Mesh
nodes :: Dict{Int64, Vector{Float64}}
node_sets :: Dict{Symbol, Set{Int64}}
elements :: Dict{Int64, Vector{Int64}}
element_types :: Dict{Int64, Symbol}
element_codes :: Dict{Int64, Symbol}
element_sets :: Dict{Symbol, Set{Int64}}
surface_sets :: Dict{Symbol, Vector{Tuple{Int64, Symbol}}}
nodes :: Dict{Int, Vector{Float64}}
node_sets :: Dict{Symbol, Set{Int}}
elements :: Dict{Int, Vector{Int}}
element_types :: Dict{Int, Symbol}
element_codes :: Dict{Int, Symbol}
element_sets :: Dict{Symbol, Set{Int}}
surface_sets :: Dict{Symbol, Vector{Tuple{Int, Symbol}}}
surface_types :: Dict{Symbol, Symbol}
end
@@ -35,7 +35,7 @@ function add_node!(mesh::Mesh, nid::Int, ncoords::Vector{Float64})
mesh.nodes[nid] = ncoords
end
function add_nodes!(mesh::Mesh, nodes::Dict{Int64, Vector{Float64}})
function add_nodes!(mesh::Mesh, nodes::Dict{Int, Vector{Float64}})
for (nid, ncoords) in nodes
add_node!(mesh, nid, ncoords)
end
@@ -43,17 +43,28 @@ end
function add_node_to_node_set!(mesh::Mesh, set_name, nids...)
if !haskey(mesh.node_sets, set_name)
mesh.node_sets[set_name] = Set{Int64}()
mesh.node_sets[set_name] = Set{Int}()
end
push!(mesh.node_sets[set_name], nids...)
return
end
function add_element!(mesh::Mesh, elid::Int, eltype::Symbol, connectivity::Vector{Int64})
""" Create a new node set from nodes in element set. """
function create_node_set_from_element_set!(mesh::Mesh, set_name)
node_ids = Set{Int}()
for elid in mesh.element_sets[set_name]
push!(node_ids, mesh.elements[elid]...)
end
mesh.node_sets[set_name] = node_ids
return
end
function add_element!(mesh::Mesh, elid::Int, eltype::Symbol, connectivity::Vector{Int})
mesh.elements[elid] = connectivity
mesh.element_types[elid] = eltype
end
function add_elements!(mesh::Mesh, elements::Dict{Int64, Tuple{Symbol, Vector{Int64}}})
function add_elements!(mesh::Mesh, elements::Dict{Int, Tuple{Symbol, Vector{Int}}})
for (elid, (eltype, elcon)) in elements
add_element!(mesh, elid, eltype, elcon)
end
@@ -61,7 +72,7 @@ end
function add_element_to_element_set!(mesh::Mesh, set_name, elids...)
if !haskey(mesh.element_sets, set_name)
mesh.element_sets[set_name] = Set{Int64}()
mesh.element_sets[set_name] = Set{Int}()
end
push!(mesh.element_sets[set_name], elids...)
end
@@ -76,7 +87,7 @@ function copy(mesh::Mesh)
return mesh2
end
function filter_by_element_id(mesh::Mesh, element_ids::Vector{Int64})
function filter_by_element_id(mesh::Mesh, element_ids::Vector{Int})
mesh2 = copy(mesh)
mesh2.elements = Dict()
for elid in element_ids
@@ -113,7 +124,7 @@ function create_elements(mesh::Mesh, element_sets::Symbol...; element_type=nothi
if isempty(element_sets)
element_ids = collect(keys(mesh.elements))
else
element_ids = Set{Int64}()
element_ids = Set{Int}()
for set_name in element_sets
element_ids = union(element_ids, mesh.element_sets[set_name])
end
@@ -135,7 +146,7 @@ end
""" find npts nearest nodes from mesh and return id numbers as list. """
function find_nearest_nodes(mesh::Mesh, coords::Vector, npts=1)
dist = Dict{Int64, Float64}()
dist = Dict{Int, Float64}()
for (nid, c) in mesh.nodes
dist[nid] = norm(coords-c)
end
@@ -166,7 +177,7 @@ function reorder_element_connectivity!(mesh::Mesh, mapping::Dict{Symbol, Vector{
end
end
function JuliaFEM.Problem{P<:FieldProblem}(mesh::Mesh, ::Type{P}, name::AbstractString, dimension::Int64)
function JuliaFEM.Problem{P<:FieldProblem}(mesh::Mesh, ::Type{P}, name::AbstractString, dimension::Int)
problem = Problem(P, name, dimension)
problem.elements = create_elements(mesh, name)
return problem
@@ -177,18 +188,3 @@ function JuliaFEM.Problem{P<:BoundaryProblem}(mesh::Mesh, ::Type{P}, name, dimen
problem.elements = create_elements(mesh, name)
return problem
end
"""
Swap surface element connectivity s.t. normals point outward
"""
function check_orientation!
# TODO
end
"""
Partition model using METIS
"""
function partition_model!
# TODO
end
src/problems_mortar.jl
+194 -3
View File
@@ -42,12 +42,14 @@ type Mortar <: BoundaryProblem
linear_surface_elements :: Bool
split_quadratic_slave_elements :: Bool
split_quadratic_master_elements :: Bool
alpha :: Float64
drop_tolerance :: Float64
store_fields :: Vector{Symbol}
end
function Mortar()
default_fields = []
return Mortar(-1, false, false, false, false, Inf, true, true, true, default_fields)
return Mortar(-1, false, false, false, false, Inf, true, true, true, 0.0, 1.0e-9, default_fields)
end
function get_unknown_field_name(problem::Problem{Mortar})
@@ -63,13 +65,202 @@ function get_formulation_type(problem::Problem{Mortar})
end
function assemble!(problem::Problem{Mortar}, time::Float64)
if length(problem.elements) == 0
warn("No elements defined in interface $(problem.name), this will result empty assembly!")
return
end
if problem.properties.dimension == -1
problem.properties.dimension = dim = size(first(problem.elements), 1)
info("assuming dimension of mesh tie surface is $dim")
info("if this is wrong set is manually using problem.properties.dimension")
info("Assuming dimension of mesh tie surface is $dim. If this is wrong set is manually using problem.properties.dimension")
end
dimension = Val{problem.properties.dimension}
use_forwarddiff = Val{problem.properties.use_forwarddiff}
assemble!(problem, time, dimension, use_forwarddiff)
end
""" Given a CCW ordered set of vertices, calculate area of polygon.
Examples
--------
julia> P = Vector[[1/3, 5/12, 1/2], [1/3, 1/2, 1/2], [1/2, 1/2, 1/2], [1/2, 1/3, 1/2], [5/12, 1/3, 1/2]]
5-element Array{Array{T,1},1}:
[0.333333,0.416667,0.5]
[0.333333,0.5,0.5]
[0.5,0.5,0.5]
[0.5,0.333333,0.5]
[0.416667,0.333333,0.5]
julia> A = calculate_polygon_area(P)
0.02430555555555556
julia> isapprox(A, 7/288)
true
"""
function calculate_polygon_area(P)
N_P = length(P)
A = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P])
return A
end
""" Function to print useful debug information from interface to find bugs. """
function diagnose_interface(problem::Problem{Mortar}, time::Float64)
info("Diagnosing Mortar interface...")
props = problem.properties
field_dim = get_unknown_field_dimension(problem)
field_name = get_parent_field_name(problem)
slave_elements = get_slave_elements(problem)
I_area = 0.0
if props.split_quadratic_slave_elements
info("props.split_quadratic_slave_elements = true")
if !props.linear_surface_elements
warn("Mortar3D: split_quadratic_surfaces = true and linear_surface_elements = false maybe have unexpected behavior")
end
slave_elements = split_quadratic_elements(slave_elements, time)
end
info("Number of slave elements in interface: $(length(slave_elements))")
# 1. calculate nodal normals and tangents for slave element nodes j ∈ S
normals = calculate_normals(slave_elements, time, Val{2};
rotate_normals=props.rotate_normals)
update!(slave_elements, "normal", time => normals)
S_areas = []
C_areas = []
P_areas = []
for slave_element in slave_elements
info(repeat("-", 80))
info("Processing slave element $(slave_element.id), type = $(get_element_type(slave_element))")
info(repeat("-", 80))
S_area = 0.0
S_area_in_contact = 0.0
for ip in get_integration_points(slave_element)
S_area += ip.weight*slave_element(ip, time, Val{:detJ})
end
info("Total area of slave element = $S_area")
if props.linear_surface_elements
info("Converting slave element to linear surface element")
slave_element = convert_to_linear_element(slave_element)
end
slave_element_nodes = get_connectivity(slave_element)
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])
# project slave nodes to auxiliary plane (x0, Q)
xi = mean(get_reference_coordinates(slave_element))
N = vec(get_basis(slave_element, xi, time))
x0 = N*X1
n0 = 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)
info("Slave element $(slave_element.id) vertices in auxiliary plane: $S")
# 3. loop all master elements
master_elements = slave_element("master elements", time)
if props.split_quadratic_master_elements
master_elements = split_quadratic_elements(master_elements, time)
end
for master_element in master_elements
if props.linear_surface_elements
master_element = convert_to_linear_element(master_element)
end
master_element_nodes = get_connectivity(master_element)
nm = length(master_element)
X2 = master_element("geometry", time)
if norm(mean(X1) - mean(X2)) > problem.properties.distval
# elements are "far enough"
continue
end
# 3.1 project master nodes to auxiliary plane and create polygon clipping
M = Vector[project_vertex_to_auxiliary_plane(X2[i], x0, n0) for i=1:nm]
check_orientation!(M, n0)
P = get_polygon_clip(S, M, n0)
if length(P) < 3
if length(P) == 0
continue
end
if length(P) == 1
info("length(P) == 1, shared vertex")
end
if length(P) == 2
info("length(P) == 2, shared edge")
end
continue
end
info("Master element $(master_element.id) vertices in auxiliary plane = $M")
check_orientation!(P, n0)
P_area_ = calculate_polygon_area(P)
info("Polygon clip found, P=$P, N_P = $(length(P)), area of polygon = $P_area_")
if isapprox(P_area_, 0.0)
error("Polygon P has zero area: $P_area_")
end
P_area = 0.0
C0 = calculate_centroid(P)
info("Centroid of polygon = $C0")
# 4. loop integration cells
all_cells = get_cells(P, C0)
info("Polygon is splitted to $(length(all_cells)) integration cells.")
for (cell_id, cell) in enumerate(all_cells)
C_area = 0.0
virtual_element = Element(Tri3, Int[])
update!(virtual_element, "geometry", cell)
# 5. loop integration point of integration cell
for ip in get_integration_points(virtual_element, 3)
N = vec(get_basis(virtual_element, ip, time))
detJ = virtual_element(ip, time, Val{:detJ})
w = ip.weight*detJ
# project gauss point from auxiliary plane to master and slave element
x_gauss = virtual_element("geometry", ip, time)
xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time)
xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, X2, time)
C_area += w
end # integration points done
info("Cell $cell_id has area of $C_area")
P_area += C_area
push!(C_areas, C_area)
end # integration cells done
if !isapprox(P_area, P_area_)
error("P_area = $P_area, should be $P_area_")
end
S_area_in_contact += P_area
push!(P_areas, P_area)
end # master elements done
S_perc = S_area_in_contact / S_area * 100.0
push!(S_areas, S_area_in_contact)
info("Area of slave element in contact: $S_area_in_contact, it's $S_perc % of total element area")
I_area += S_area_in_contact
end # slave elements done, contact virtual work ready
info("Area of interface: $I_area")
info("Smallest cell area: $(minimum(C_areas))")
info("Smallest polygon area: $(minimum(P_areas))")
info("Smallest slave element area in contact: $(minimum(S_areas))")
end
src/problems_mortar_3d.jl
+524 -190
View File
@@ -112,7 +112,7 @@ function get_polygon_clip(xs, xm, n)
# 1. test is master point inside slave, if yes, add to clip
for i=1:nm
if vertex_inside_polygon(xm[i], xs)
debug("1. $(xm[i]) inside S -> push")
# debug("1. $(xm[i]) inside S -> push")
push!(P, xm[i])
end
end
@@ -121,7 +121,7 @@ function get_polygon_clip(xs, xm, n)
for i=1:ns
if vertex_inside_polygon(xs[i], xm)
approx_in(xs[i], P) && continue
debug("2. $(xs[i]) inside M -> push")
# debug("2. $(xs[i]) inside M -> push")
push!(P, xs[i])
end
end
@@ -144,7 +144,7 @@ function get_polygon_clip(xs, xm, n)
#info("t=$t, q=$q, q ∈ xm ? $(vertex_inside_polygon(q, xm))")
if vertex_inside_polygon(q, xm)
approx_in(q, P) && continue
debug("3. $q inside M -> push")
# debug("3. $q inside M -> push")
push!(P, q)
end
end
@@ -228,12 +228,39 @@ function calculate_normals(elements, time, ::Type{Val{2}}; rotate_normals=false)
return normals
end
""" Given polygon P and normal direction n, check that polygon vertices are
ordered in counter clock wise direction with respect to surface normal and
sort if necessary. It is assumed that polygon is convex.
Examples
--------
Unit triangle, normal in z-direction:
julia> P = Vector[[0.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 0.0, 0.0]]
3-element Array{Array{T,1},1}:
[0.0,0.0,0.0]
[0.0,1.0,0.0]
[1.0,0.0,0.0]
julia> n = [0.0, 0.0, 1.0]
3-element Array{Float64,1}:
0.0
0.0
1.0
julia> check_orientation!(P, n)
3-element Array{Array{T,1},1}:
[1.0,0.0,0.0]
[0.0,0.0,0.0]
[0.0,1.0,0.0]
"""
function check_orientation!(P, n)
C = mean(P)
np = length(P)
s = [dot(n, cross(P[i]-C, P[mod(i+1,np)+1]-C)) for i=1:np]
all(s .< 0) && return
debug("polygon not in ccw order, fixing")
# debug("polygon not in ccw order, fixing")
# project points to new orthogonal basis Q and sort there
t1 = (P[1]-C)/norm(P[1]-C)
t2 = cross(n, t1)
@@ -274,10 +301,9 @@ function split_quadratic_element(element::Element{Tri6}, time::Float64)
u = element("displacement", time)
update!(new_element, "displacement", time => u[elmap])
end
#n = element("normal", time)
#update!(new_element, "normal", time => n[elmap])
if haskey(element, "master elements")
update!(new_element, "master elements", time => element("master elements", time))
if haskey(element, "normal")
n = element("normal", time)
update!(new_element, "normal", time => n[elmap])
end
push!(new_elements, new_element)
end
@@ -298,70 +324,62 @@ function split_quadratic_elements(elements::Vector, time::Float64)
end
n1 = length(elements)
n2 = length(new_elements)
info("Splitted $n1 (maybe quadratic) elements to $n2 (linear) sub-elements")
if n1 != n2
info("Splitted $n1 elements to $n2 (linear) sub-elements")
end
return new_elements
end
function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{Val{false}})
""" Assemble linear surface element to problem.
Dual basis is constructed such that partially integrated slave segments are taken into account in a proper way.
Notes
-----
For full integrated slave element, coefficient matrix for Tri3 is
Ae = [3.0 -1.0 -1.0; -1.0 3.0 -1.0; -1.0 -1.0 3.0]
References
----------
[Popp2013] Popp, Alexander, et al. "Improved robustness and consistency of 3D contact algorithms based on a dual mortar approach." Computer Methods in Applied Mechanics and Engineering 264 (2013): 67-80.
"""
function assemble!{E<:Union{Tri3, Quad4}}(problem::Problem{Mortar}, slave_element::Element{E}, time::Real; first_slave_element=false)
props = problem.properties
field_dim = get_unknown_field_dimension(problem)
field_name = get_parent_field_name(problem)
slave_elements = get_slave_elements(problem)
area = 0.0
if props.split_quadratic_slave_elements
if !props.linear_surface_elements
warn("Mortar3D: split_quadratic_surfaces = true and linear_surface_elements = false maybe have unexpected behavior")
end
slave_elements = split_quadratic_elements(slave_elements, time)
end
slave_element_nodes = get_connectivity(slave_element)
nsl = length(slave_element)
X1 = slave_element("geometry", time)
n1 = slave_element("normal", time)
# 1. calculate nodal normals and tangents for slave element nodes j ∈ S
normals = calculate_normals(slave_elements, time, Val{2};
rotate_normals=props.rotate_normals)
update!(slave_elements, "normal", time => normals)
# project slave nodes to auxiliary plane (x0, Q)
xi = mean(get_reference_coordinates(slave_element))
first_slave_element && debug("midpoint xi = $xi")
N = vec(get_basis(slave_element, xi, time))
x0 = N*X1
n0 = N*n1
S = Vector[project_vertex_to_auxiliary_plane(X1[i], x0, n0) for i=1:nsl]
# 2. loop all slave elements
first_slave_element = true
master_elements = slave_element("master elements", time)
for slave_element in slave_elements
if props.dual_basis
if props.linear_surface_elements
slave_element = convert_to_linear_element(slave_element)
end
slave_element_nodes = get_connectivity(slave_element)
nsl = length(slave_element)
X1 = slave_element("geometry", time)
n1 = Field([normals[j] for j in slave_element_nodes])
# project slave nodes to auxiliary plane (x0, Q)
xi = mean(get_reference_coordinates(slave_element))
first_slave_element && debug("midpoint xi = $xi")
N = vec(get_basis(slave_element, xi, time))
x0 = N*X1
n0 = N*n1
S = Vector[project_vertex_to_auxiliary_plane(X1[i], x0, n0) for i=1:nsl]
# 3. loop all master elements
master_elements = slave_element("master elements", time)
if props.split_quadratic_master_elements
master_elements = split_quadratic_elements(master_elements, time)
end
debug("Creating dual basis for element $(slave_element.id)")
De = zeros(nsl, nsl)
Me = zeros(nsl, nsl)
for master_element in master_elements
if props.linear_surface_elements
master_element = convert_to_linear_element(master_element)
end
master_element_nodes = get_connectivity(master_element)
nm = length(master_element)
X2 = master_element("geometry", time)
if norm(mean(X1) - mean(X2)) > problem.properties.distval
# elements are "far enough"
continue
end
@@ -373,170 +391,486 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{
N_P = length(P)
P_area = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P])
if first_slave_element
debug("Polygon clip info for first slave element:")
debug("S = $S")
debug("M = $M")
debug("P = $P")
debug("N_P = $N_P")
debug("P_area = $P_area")
end
if isapprox(P_area, 0.0)
info("Polygon P has zero area: $P_area")
continue
end
C0 = calculate_centroid(P)
#=
if isnan(C0[1])
info("C0 = $C0")
info("P = $P")
info("S = $S")
info("M = $M")
info("n0 = $n0")
error("Calculation of centroid of polygon clip P failed.")
end
=#
De = zeros(nsl, nsl)
Me = zeros(nsl, nm)
ge = zeros(field_dim*nsl)
# 4. loop integration cells
C0 = calculate_centroid(P)
all_cells = get_cells(P, C0)
for cell in all_cells
virtual_element = Element(Tri3, Int[])
update!(virtual_element, "geometry", cell)
#x_cell = Field(cell)
# construct bi-orthogonal basis
nnodes = length(slave_element)
if props.dual_basis
De = zeros(nnodes, nnodes)
Me = zeros(nnodes, nnodes)
for ip in get_integration_points(virtual_element, 3)
x_gauss = nothing
#try
x_gauss = virtual_element("geometry", ip, time)
xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time)
detJ = virtual_element(ip, time, Val{:detJ})
w = ip.weight*detJ
N1 = vec(get_basis(slave_element, xi_s, time))
De += w*diagm(vec(N1))
Me += w*N1*N1'
#catch
# info("Failed to construct bi-orthogonal basis: cannot project vertex from auxiliary plane back to sufface.")
# info("x_gauss = $x_gauss")
# info("cell = $cell")
# info("C0 = $C0")
# info("P = $P")
# info("S = $S")
# info("M = $M")
# info("n0 = $n0")
# rethrow()
#end
end
Ae = De*inv(Me)
else
Ae = eye(nnodes)
end
# 5. loop integration point of integration cell
for ip in get_integration_points(virtual_element, 3)
N = vec(get_basis(virtual_element, ip, time))
#dN = vec(get_dbasis(virtual_element, ip, time))
#JC = transpose(sum([kron(dNC[:,j], x_cell[j]') for j=1:length(x_cell)]))
#wC = ip.weight*norm(cross(JC[:,1], JC[:,2]))
detJ = virtual_element(ip, time, Val{:detJ})
w = ip.weight*detJ
# project gauss point from auxiliary plane to master and slave element
#x_gauss = N*x_cell
x_gauss = virtual_element("geometry", ip, time)
#=
if isnan(x_gauss[1])
info("is nan")
info("x_gauss = $x_gauss")
info("cell = $cell")
info("C0 = $C0")
info("P = $P")
info("S = $S")
info("M = $M")
info("n0 = $n0")
error("nan, unable to continue")
end
=#
xi_s = nothing
xi_m = nothing
alpha = nothing
#try
xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time)
xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, X2, time)
#catch
# info("projecting vertex back to surface has failed.")
# info("x_gauss = $x_gauss")
# info("cell = $cell")
# info("C0 = $C0")
# info("P = $P")
# info("S = $S")
# info("M = $M")
# info("n0 = $n0")
# rethrow()
#end
# add contributions
N1 = vec(get_basis(slave_element, xi_s, time))
N2 = vec(get_basis(master_element, xi_m, time))
Phi = Ae*N1
De += w*Phi*N1'
Me += w*Phi*N2'
if props.adjust
u1 = slave_element("displacement", time)
u2 = master_element("displacement", time)
x_s = N1*(X1+u1)
x_m = N2*(X2+u2)
ge += w*vec((x_m-x_s)*Phi')
end
area += w
end # integration points done
xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time)
N1 = slave_element(xi_s, time)
De += w*diagm(vec(N1))
Me += w*N1'*N1
end
end # integration cells done
# 6. add contribution to contact virtual work
sdofs = get_gdofs(problem, slave_element)
mdofs = get_gdofs(problem, master_element)
for i=1:field_dim
lsdofs = sdofs[i:field_dim:end]
lmdofs = mdofs[i:field_dim:end]
add!(problem.assembly.C1, lsdofs, lsdofs, De)
add!(problem.assembly.C1, lsdofs, lmdofs, -Me)
add!(problem.assembly.C2, lsdofs, lsdofs, De)
add!(problem.assembly.C2, lsdofs, lmdofs, -Me)
end # master elements done
Ae = De*inv(Me)
info("Dual basis coefficient matrix: $Ae")
else
Ae = eye(nsl)
end
for master_element in master_elements
master_element_nodes = get_connectivity(master_element)
nm = length(master_element)
X2 = master_element("geometry", time)
if norm(mean(X1) - mean(X2)) > problem.properties.distval
continue
end
# 3.1 project master nodes to auxiliary plane and create polygon clipping
M = Vector[project_vertex_to_auxiliary_plane(X2[i], x0, n0) for i=1:nm]
P = get_polygon_clip(S, M, n0)
length(P) < 3 && continue # no clipping or shared edge (no volume)
check_orientation!(P, n0)
N_P = length(P)
P_area = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P])
if first_slave_element
debug("Polygon clip info for first slave element:")
debug("S = $S")
debug("M = $M")
debug("P = $P")
debug("N_P = $N_P")
debug("P_area = $P_area")
end
if isapprox(P_area, 0.0)
info("Polygon P has zero area: $P_area")
continue
end
C0 = calculate_centroid(P)
De = zeros(nsl, nsl)
Me = zeros(nsl, nm)
ge = zeros(field_dim*nsl)
# 4. loop integration cells
all_cells = get_cells(P, C0)
for cell in all_cells
virtual_element = Element(Tri3, Int[])
update!(virtual_element, "geometry", cell)
# 5. loop integration point of integration cell
for ip in get_integration_points(virtual_element, 3)
detJ = virtual_element(ip, time, Val{:detJ})
w = ip.weight*detJ
# project gauss point from auxiliary plane to master and slave element
x_gauss = virtual_element("geometry", ip, time)
xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time)
xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, X2, time)
# add contributions
N1 = vec(get_basis(slave_element, xi_s, time))
N2 = vec(get_basis(master_element, xi_m, time))
Phi = Ae*N1
# Phi = [3.0-4.0*xi_s[1]-4.0*xi_s[2], 4.0*xi_s[1]-1.0, 4.0*xi_s[2]-1.0]
De += w*Phi*N1'
Me += w*Phi*N2'
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)
ge += w*vec((x_m-x_s)*Phi')
end
area += w
end # integration points done
end # integration cells done
# 6. add contribution to contact virtual work
sdofs = get_gdofs(problem, slave_element)
mdofs = get_gdofs(problem, master_element)
for i=1:field_dim
lsdofs = sdofs[i:field_dim:end]
lmdofs = mdofs[i:field_dim:end]
add!(problem.assembly.C1, lsdofs, lsdofs, De)
add!(problem.assembly.C1, lsdofs, lmdofs, -Me)
add!(problem.assembly.C2, lsdofs, lsdofs, De)
add!(problem.assembly.C2, lsdofs, lmdofs, -Me)
end
add!(problem.assembly.g, sdofs, ge)
end # master elements done
return area
end
""" Assemble quadratic surface element to problem.
In polygon clipping element is divided to linear sub-elements proposed in [Puso2008].
References
----------
[Puso2008] Puso, Michael A., T. A. Laursen, and Jerome Solberg. "A segment-to-segment mortar contact method for quadratic elements and large deformations." Computer Methods in Applied Mechanics and Engineering 197.6 (2008): 555-566.
[Popp1012] Popp, Alexander, et al. "Dual quadratic mortar finite element methods for 3D finite deformation contact." SIAM Journal on Scientific Computing 34.4 (2012): B421-B446.
"""
function assemble!{E<:Union{Tri6}}(problem::Problem{Mortar}, slave_element::Element{E}, time::Real; first_slave_element=false)
props = problem.properties
field_dim = get_unknown_field_dimension(problem)
field_name = get_parent_field_name(problem)
area = 0.0
Xs = slave_element("geometry", time)
alp = props.alpha
if alp != 0.0
T = [
1.0 0.0 0.0 0.0 0.0 0.0
0.0 1.0 0.0 0.0 0.0 0.0
0.0 0.0 1.0 0.0 0.0 0.0
alp alp 0.0 1.0-2*alp 0.0 0.0
0.0 alp alp 0.0 1.0-2*alp 0.0
alp 0.0 alp 0.0 0.0 1.0-2*alp
]
else
T = eye(6)
end
#=
invT = [
1.0 0.0 0.0 0.0 0.0 0.0
0.0 1.0 0.0 0.0 0.0 0.0
0.0 0.0 1.0 0.0 0.0 0.0
-alp/(1-2*alp) -alp/(1-2*alp) 0.0 1/(1-2*alp) 0.0 0.0
0.0 -alp/(1-2*alp) -alp/(1-2*alp) 0.0 1/(1-2*alp) 0.0
-alp/(1-2*alp) 0.0 -alp/(1-2*alp) 0.0 0.0 1/(1-2*alp)
]
=#
if props.dual_basis
# info("Creating dual basis for element $(slave_element.id)")
nsl = length(slave_element)
De = zeros(nsl, nsl)
Me = zeros(nsl, nsl)
# split slave element to linear sub-elements and loop
for sub_slave_element in split_quadratic_element(slave_element, time)
slave_element_nodes = get_connectivity(sub_slave_element)
nsl = length(sub_slave_element)
X1 = sub_slave_element("geometry", time)
n1 = sub_slave_element("normal", time)
# create auxiliary plane
xi = mean(get_reference_coordinates(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
# project slave nodes to auxiliary plane
S = Vector[project_vertex_to_auxiliary_plane(X1[i], x0, n0) for i=1:nsl]
# 3. loop all master elements
master_elements = slave_element("master elements", time)
for master_element in master_elements
Xm = master_element("geometry", time)
if norm(mean(Xs) - mean(Xm)) > problem.properties.distval
continue
end
# split master element to linear sub-elements and loop
for sub_master_element in split_quadratic_element(master_element, time)
master_element_nodes = get_connectivity(sub_master_element)
nm = length(sub_master_element)
X2 = sub_master_element("geometry", time)
# 3.1 project master nodes to auxiliary plane
M = Vector[project_vertex_to_auxiliary_plane(X2[i], x0, n0) for i=1:nm]
# create polygon clipping P
P = get_polygon_clip(S, M, n0)
length(P) < 3 && continue # no clipping or shared edge (no volume)
check_orientation!(P, n0)
N_P = length(P)
P_area = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P])
C0 = calculate_centroid(P)
# 4. loop integration cells
all_cells = get_cells(P, C0)
for cell in all_cells
virtual_element = Element(Tri3, Int[])
update!(virtual_element, "geometry", 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)
detJ = virtual_element(ip, time, Val{:detJ})
w = ip.weight*detJ
N1 = vec(slave_element(xi_s, time)*T)
De += w*diagm(N1)
Me += w*N1*N1'
end
end # integration cells done
end # sub aster elements done
end # master elements done
end # sub slave elements done
Ae = De*inv(Me)
# info("Dual basis construction finished.")
# info("Slave element geometry = $Xs")
# info("De = $De")
# info("Me = $Me")
# info("Dual basis coefficient matrix: $Ae")
else
nsl = length(slave_element)
Ae = eye(nsl)
end
# split slave element to linear sub-elements and loop
for sub_slave_element in split_quadratic_element(slave_element, time)
slave_element_nodes = get_connectivity(sub_slave_element)
nsl = length(sub_slave_element)
X1 = sub_slave_element("geometry", time)
n1 = sub_slave_element("normal", time)
# create auxiliary plane
xi = mean(get_reference_coordinates(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
# project slave nodes to auxiliary plane
S = Vector[project_vertex_to_auxiliary_plane(X1[i], x0, n0) for i=1:nsl]
# 3. loop all master elements
master_elements = slave_element("master elements", time)
for master_element in master_elements
Xm = master_element("geometry", time)
if norm(mean(Xs) - mean(Xm)) > problem.properties.distval
continue
end
add!(problem.assembly.g, sdofs, ge)
# split master element to linear sub-elements and loop
for sub_master_element in split_quadratic_element(master_element, time)
master_element_nodes = get_connectivity(sub_master_element)
nm = length(sub_master_element)
X2 = sub_master_element("geometry", time)
# 3.1 project master nodes to auxiliary plane
M = Vector[project_vertex_to_auxiliary_plane(X2[i], x0, n0) for i=1:nm]
# create polygon clipping P
P = get_polygon_clip(S, M, n0)
length(P) < 3 && continue # no clipping or shared edge (no volume)
check_orientation!(P, n0)
N_P = length(P)
P_area = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P])
if first_slave_element
debug("Polygon clip info for first slave element:")
debug("S = $S")
debug("M = $M")
debug("P = $P")
debug("N_P = $N_P")
debug("P_area = $P_area")
end
if isapprox(P_area, 0.0)
warn("Polygon P has zero area: $P_area")
continue
end
C0 = calculate_centroid(P)
# while our polygon clipping algorithm is working in linear sub elements
# contributions is calculated using quadratic shape functions
De = zeros(length(slave_element), length(slave_element))
Me = zeros(length(slave_element), length(master_element))
ge = zeros(field_dim*length(slave_element))
# 4. loop integration cells
all_cells = get_cells(P, C0)
for cell in all_cells
virtual_element = Element(Tri3, Int[])
update!(virtual_element, "geometry", cell)
# 5. loop integration point of integration 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)
xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, Xm, time)
# add contributions
N1 = vec(slave_element(xi_s, time)*T)
N2 = vec(master_element(xi_m, time))
Phi = Ae*N1
detJ = virtual_element(ip, time, Val{:detJ})
w = ip.weight*detJ
De += w*Phi*N1'
Me += w*Phi*N2'
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)
ge += w*vec((xm-xs)*Phi')
end
area += w
end # integration points done
end # integration cells done
# 6. add contribution to contact virtual work
sdofs = get_gdofs(problem, slave_element)
mdofs = get_gdofs(problem, master_element)
for i=1:field_dim
lsdofs = sdofs[i:field_dim:end]
lmdofs = mdofs[i:field_dim:end]
add!(problem.assembly.C1, lsdofs, lsdofs, De)
add!(problem.assembly.C1, lsdofs, lmdofs, -Me)
add!(problem.assembly.C2, lsdofs, lsdofs, De)
add!(problem.assembly.C2, lsdofs, lmdofs, -Me)
end
add!(problem.assembly.g, sdofs, ge)
end # sub aster elements done
end # master elements done
end # sub slave elements done
return area
end
function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{Val{false}})
props = problem.properties
field_dim = get_unknown_field_dimension(problem)
field_name = get_parent_field_name(problem)
slave_elements = get_slave_elements(problem)
area = 0.0
#=
if props.split_quadratic_slave_elements
if !props.linear_surface_elements
warn("Mortar3D: split_quadratic_surfaces = true and linear_surface_elements = false maybe have unexpected behavior")
end
slave_elements = split_quadratic_elements(slave_elements, time)
end
=#
# 1. calculate nodal normals and tangents for slave element nodes j ∈ S
normals = calculate_normals(slave_elements, time, Val{2};
rotate_normals=props.rotate_normals)
update!(slave_elements, "normal", time => normals)
# 2. loop all slave elements
first_slave_element = true
for slave_element in slave_elements
area += assemble!(problem, slave_element, time; first_slave_element=first_slave_element)
first_slave_element = false
end # slave elements done, contact virtual work ready
if problem.properties.dual_basis
tol = 1.0e-9
debug("Dual basis is used, dropping small values for C1 & C2, tol = $tol")
C1 = sparse(problem.assembly.C1)
C2 = sparse(problem.assembly.C2)
SparseArrays.droptol!(C1, tol)
SparseArrays.droptol!(C2, tol)
problem.assembly.C1 = C1
problem.assembly.C2 = C2
C1 = sparse(problem.assembly.C1)
C2 = sparse(problem.assembly.C2)
maxdim = maximum(size(C1))
if problem.properties.alpha != 0.0
debug("mortar_3d: size C1 = ", size(C1), " max dim = $maxdim")
debug("alpha != 0.0, applying transformation D = Dh*T^-1")
alp = problem.properties.alpha
Te = [
1.0 0.0 0.0 0.0 0.0 0.0
0.0 1.0 0.0 0.0 0.0 0.0
0.0 0.0 1.0 0.0 0.0 0.0
alp alp 0.0 1.0-2*alp 0.0 0.0
0.0 alp alp 0.0 1.0-2*alp 0.0
alp 0.0 alp 0.0 0.0 1.0-2*alp
]
invTe = [
1.0 0.0 0.0 0.0 0.0 0.0
0.0 1.0 0.0 0.0 0.0 0.0
0.0 0.0 1.0 0.0 0.0 0.0
-alp/(1-2*alp) -alp/(1-2*alp) 0.0 1/(1-2*alp) 0.0 0.0
0.0 -alp/(1-2*alp) -alp/(1-2*alp) 0.0 1/(1-2*alp) 0.0
-alp/(1-2*alp) 0.0 -alp/(1-2*alp) 0.0 0.0 1/(1-2*alp)
]
# construct global transformation matrices T and invT
T = SparseMatrixCOO()
invT = SparseMatrixCOO()
for element in slave_elements
dofs = get_gdofs(problem, element)
for i=1:field_dim
ldofs = dofs[i:field_dim:end]
add!(T, ldofs, ldofs, Te)
add!(invT, ldofs, ldofs, invTe)
end
end
T = sparse(T, maxdim, maxdim, (a, b) -> b)
invT = sparse(invT, maxdim, maxdim, (a, b) -> b)
# fill diagonal
d = ones(size(T, 1))
d[get_nonzero_rows(T)] = 0.0
T += spdiagm(d)
invT += spdiagm(d)
#invT2 = sparse(inv(full(T)))
#info("invT == invT2? ", invT == invT2)
#maxabsdiff = maximum(abs(invT - invT2))
#info("max diff = $maxabsdiff")
C1 = C1*invT
C2 = C2*invT
end
tol = problem.properties.drop_tolerance
debug("Dropping small values from C1 & C2, tolerace = $tol")
SparseArrays.droptol!(C1, tol)
SparseArrays.droptol!(C2, tol)
problem.assembly.C1 = C1
problem.assembly.C2 = C2
debug("area of interface: $area")
end
src/solvers.jl
+5
View File
@@ -109,11 +109,13 @@ end
function check_for_overconstrained_dofs(solver::Solver)
overdetermined = false
constrained_dofs = Set{Int}()
all_overconstrained_dofs = Set{Int}()
boundary_problems = get_boundary_problems(solver)
for problem in boundary_problems
new_constraints = Set(problem.assembly.C2.I)
new_constraints = setdiff(new_constraints, problem.assembly.removed_dofs)
overconstrained_dofs = intersect(constrained_dofs, new_constraints)
all_overconstrained_dofs = union(all_overconstrained_dofs, overconstrained_dofs)
if length(overconstrained_dofs) != 0
warn("problem is overconstrained, finding overconstrained dofs... ")
overdetermined = true
@@ -133,6 +135,8 @@ function check_for_overconstrained_dofs(solver::Solver)
constrained_dofs = union(constrained_dofs, new_constraints)
end
if overdetermined
warn("List of all overconstrained dofs:")
warn(sort(collect(all_overconstrained_dofs)))
error("problem is overconstrained, not continuing to solution.")
end
return true
@@ -276,6 +280,7 @@ function solve!(solver::Solver, K, C1, C2, D, f, g, u, la, ::Type{Val{3}})
u[:] = x[1:solver.ndofs]
la[:] = x[solver.ndofs+1:end]
return true
end