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

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 5631f2895b
commit 4a471b5cb7
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
test/test_mortar_3d_polygon_clip.jl
+15 -14
View File
@@ -1,10 +1,11 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using JuliaFEM
using JuliaFEM: get_polygon_clip, calculate_polygon_area
using JuliaFEM.Testing
@testset "polygon clip case 1" begin
@testset "polygon clipping" begin
S = Vector[
[0.375, 0.0, 0.5],
[0.6, 0.0, 0.5],
@@ -15,17 +16,9 @@ using JuliaFEM.Testing
[0.375, 0.25, 0.5]]
n0 = [0.0, 0.0, 1.0]
P = get_polygon_clip(S, M, n0)
P_expected = Vector{Float64}[
[0.500, 0.0, 0.5],
[0.375, 0.0, 0.5],
[0.4375, 0.125, 0.5]]
@test length(P) == length(P_expected)
for (Pi, Pj) in zip(P, P_expected)
@test isapprox(Pi, Pj)
end
end
@testset "polygon clip case 2" begin
@test length(P) == 3
@test isapprox(calculate_polygon_area(P), 1/128)
S = Vector[
[0.25, 0.0, 0.5],
[0.75, 0.0, 0.5],
@@ -37,5 +30,13 @@ end
n0 = [0.0, 0.0, 1.0]
P = get_polygon_clip(S, M, n0)
@test length(P) == 3
end
@test isapprox(calculate_polygon_area(P), 1/48)
# visually inspected
Xs = Vector[[0.0, 0.0, 0.5], [1.0, 0.0, 0.5], [0.0, 1.0, 0.5]]
Xm = Vector[[-0.25, 0.50, 0.5], [0.50, -0.25, 0.5], [0.75,0.75, 0.5]]
P_ = Vector[[0.65,0.35,0.0], [0.5625,0.0,0.0], [0.25,0.0,0.0],
[0.0,0.25,0.0], [0.0,0.5625,0.0], [0.35,0.65,0.0]]
P = get_polygon_clip(Xs, Xm, [0.0, 0.0, 1.0])
@test length(P) == length(P_)
end
test/test_problems_mortar_3d.jl
+428 -11
View File
@@ -7,32 +7,40 @@ using JuliaFEM.Postprocess
using JuliaFEM.Testing
using JuliaFEM.Abaqus: create_surface_elements
@testset "test that interface transfers constant field without error" begin
meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/block_3d.med"
mesh = aster_read_mesh(meshfile)
### temperature patch tests, sl tet4, dl tet4, sl tet10, dl tet 10
upper = Problem(Heat, "upper", 1)
tet4_meshfile = "test_problems_mortar_3d/tet4.inp"
tet10_meshfile = "test_problems_mortar_3d/tet10.inp"
@testset "patch test temperature + abaqus inp + tet4" begin
mesh = abaqus_read_mesh(tet4_meshfile)
upper = Problem(Heat, "UPPER", 1)
upper.elements = create_elements(mesh, "UPPER")
update!(upper, "temperature thermal conductivity", 1.0)
lower = Problem(Heat, "lower", 1)
lower = Problem(Heat, "LOWER", 1)
lower.elements = create_elements(mesh, "LOWER")
update!(lower, "temperature thermal conductivity", 1.0)
bc_upper = Problem(Dirichlet, "upper boundary", 1, "temperature")
bc_upper.elements = create_elements(mesh, "UPPER_TOP")
bc_upper = Problem(Dirichlet, "UPPER_TOP", 1, "temperature")
bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper, "temperature 1", 0.0)
bc_lower = Problem(Dirichlet, "lower boundary", 1, "temperature")
bc_lower.elements = create_elements(mesh, "LOWER_BOTTOM")
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 1, "temperature")
bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower, "temperature 1", 1.0)
interface = Problem(Mortar, "interface between upper and lower block", 1, "temperature")
interface_slave_elements = create_elements(mesh, "LOWER_TOP")
interface_master_elements = create_elements(mesh, "UPPER_BOTTOM")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
interface.elements = [interface_master_elements; interface_slave_elements]
JuliaFEM.diagnose_interface(interface, 0.0)
solver = LinearSolver(upper, lower, bc_upper, bc_lower, interface)
solver.xdmf = Xdmf("sl_lin_temp_results")
solver()
node_ids, temperature = get_nodal_vector(interface.elements, "temperature", 0.0)
@@ -42,4 +50,413 @@ using JuliaFEM.Abaqus: create_surface_elements
info("minT = $minT, maxT = $maxT")
@test isapprox(minT, 0.5)
@test isapprox(maxT, 0.5)
#=
initialize!(solver)
assemble!(solver)
M, K, Kg, f, fg = get_field_assembly(solver)
Kb, C1, C2, D, fb, g = get_boundary_assembly(solver)
K = K + Kg + Kb
f = f + fg + fb
K = 1/2*(K + K')
M = 1/2*(M + M')
=#
end
@testset "patch test temperature + abaqus inp + tet4 + dual basis + adjust" begin
mesh = abaqus_read_mesh(tet4_meshfile)
upper = Problem(Heat, "UPPER", 1)
upper.elements = create_elements(mesh, "UPPER")
update!(upper, "temperature thermal conductivity", 1.0)
lower = Problem(Heat, "LOWER", 1)
lower.elements = create_elements(mesh, "LOWER")
update!(lower, "temperature thermal conductivity", 1.0)
bc_upper = Problem(Dirichlet, "UPPER_TOP", 1, "temperature")
bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper, "temperature 1", 0.0)
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 1, "temperature")
bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower, "temperature 1", 1.0)
interface = Problem(Mortar, "interface between upper and lower block", 1, "temperature")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
interface.elements = [interface_master_elements; interface_slave_elements]
interface.properties.dual_basis = true
#interface.properties.adjust = true
JuliaFEM.diagnose_interface(interface, 0.0)
solver = LinearSolver(upper, lower, bc_upper, bc_lower, interface)
solver.xdmf = Xdmf("dl_lin_temp_results")
solver()
node_ids, temperature = get_nodal_vector(interface.elements, "temperature", 0.0)
T = [t[1] for t in temperature]
minT = minimum(T)
maxT = maximum(T)
info("minT = $minT, maxT = $maxT")
@test isapprox(minT, 0.5)
@test isapprox(maxT, 0.5)
end
@testset "patch test temperature + abaqus inp + tet10, quadratic surface elements" begin
mesh = abaqus_read_mesh(tet10_meshfile)
upper = Problem(Heat, "UPPER", 1)
upper.elements = create_elements(mesh, "UPPER")
update!(upper, "temperature thermal conductivity", 1.0)
lower = Problem(Heat, "LOWER", 1)
lower.elements = create_elements(mesh, "LOWER")
update!(lower, "temperature thermal conductivity", 1.0)
bc_upper = Problem(Dirichlet, "UPPER_TOP", 1, "temperature")
bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper, "temperature 1", 0.0)
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 1, "temperature")
bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower, "temperature 1", 1.0)
interface = Problem(Mortar, "interface between upper and lower block", 1, "temperature")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
interface.elements = [interface_master_elements; interface_slave_elements]
interface.properties.linear_surface_elements = false
interface.properties.split_quadratic_slave_elements = false
interface.properties.split_quadratic_master_elements = false
interface.properties.alpha = 0.0
# JuliaFEM.diagnose_interface(interface, 0.0)
solver = LinearSolver(upper, lower, bc_upper, bc_lower, interface)
solver.xdmf = Xdmf("sl_quad_temp_results")
solver()
node_ids, temperature = get_nodal_vector(interface.elements, "temperature", 0.0)
T = [t[1] for t in temperature]
minT = minimum(T)
maxT = maximum(T)
stdT = std(T)
info("minT = $minT, maxT = $maxT, stdT = $stdT")
@test maxT - minT < 1.0e-10
@test isapprox(stdT, 0.0; atol=1.0e-10)
end
@testset "patch test temperature + abaqus inp + tet10 + quadratic surface elements + dual basis + alpha=0.2" begin
mesh = abaqus_read_mesh(tet10_meshfile)
upper = Problem(Heat, "UPPER", 1)
upper.elements = create_elements(mesh, "UPPER")
update!(upper, "temperature thermal conductivity", 1.0)
lower = Problem(Heat, "LOWER", 1)
lower.elements = create_elements(mesh, "LOWER")
update!(lower, "temperature thermal conductivity", 1.0)
bc_upper = Problem(Dirichlet, "UPPER_TOP", 1, "temperature")
bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper, "temperature 1", 0.0)
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 1, "temperature")
bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower, "temperature 1", 1.0)
interface = Problem(Mortar, "interface between upper and lower block", 1, "temperature")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
interface.elements = [interface_master_elements; interface_slave_elements]
interface.properties.linear_surface_elements = false
interface.properties.split_quadratic_slave_elements = false
interface.properties.split_quadratic_master_elements = false
interface.properties.dual_basis = true
interface.properties.alpha = 0.2
solver = LinearSolver(upper, lower, bc_upper, bc_lower, interface)
solver.xdmf = Xdmf("dl_quad_temp_results")
solver()
node_ids, temperature = get_nodal_vector(interface.elements, "temperature", 0.0)
#node_ids, temperature = get_nodal_vector(interface_slave_elements, "temperature", 0.0)
#=
for (j, (nid, T)) in enumerate(zip(node_ids, temperature))
info("$j: $nid -> $(T[1])")
j == 10 && break
end
=#
T = [t[1] for t in temperature]
minT = minimum(T)
maxT = maximum(T)
stdT = std(T)
info("minT = $minT, maxT = $maxT, stdT = $stdT")
@test maxT - minT < 1.0e-10
@test isapprox(stdT, 0.0; atol=1.0e-10)
end
### displacement patch tests, sl tet4, dl tet4, sl tet10, dl tet 10
@testset "patch test displacement + abaqus inp + tet4 + adjust" begin
mesh = abaqus_read_mesh(tet4_meshfile)
# modify mesh a bit, find all nodes in elements in element set UPPER and put 0.2 to X3 to test adjust
JuliaFEM.Preprocess.create_node_set_from_element_set!(mesh, :UPPER)
for nid in mesh.node_sets[:UPPER]
mesh.nodes[nid][3] += 0.2
end
upper = Problem(Elasticity, "UPPER", 3)
upper.elements = create_elements(mesh, "UPPER")
update!(upper, "youngs modulus", 288.0)
update!(upper, "poissons ratio", 1/3)
lower = Problem(Elasticity, "LOWER", 3)
lower.elements = create_elements(mesh, "LOWER")
update!(lower, "youngs modulus", 288.0)
update!(lower, "poissons ratio", 1/3)
bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement")
bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper, "displacement 3", 0.0)
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement")
bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower, "displacement 3", 0.0)
bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement")
bc_sym13.elements = [create_surface_elements(mesh, "LOWER_SYM13"); create_surface_elements(mesh, "UPPER_SYM13")]
update!(bc_sym13, "displacement 2", 0.0)
bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement")
bc_sym23.elements = [create_surface_elements(mesh, "LOWER_SYM23"); create_surface_elements(mesh, "UPPER_SYM23")]
update!(bc_sym23, "displacement 1", 0.0)
interface = Problem(Mortar, "LOWER_TO_UPPER", 3, "displacement")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
interface.elements = [interface_slave_elements; interface_master_elements]
append!(interface.assembly.removed_dofs, [1316, 1319, 1358, 1387, 1388, 1492, 1597, 1627])
interface.properties.linear_surface_elements = false
interface.properties.split_quadratic_slave_elements = false
interface.properties.split_quadratic_master_elements = false
interface.properties.adjust = true
interface.properties.dual_basis = false
solver = LinearSolver(upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface)
solver.xdmf = Xdmf("sl_lin_disp_results")
solver()
node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0)
node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0)
u3 = [u[3] for u in displacement]
maxabsu3 = maximum(abs(u3))
stdabsu3 = std(abs(u3))
info("tet10 block: max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-10)
end
@testset "patch test displacement + abaqus inp + tet4 + adjust + dual basis" begin
mesh = abaqus_read_mesh(tet4_meshfile)
# modify mesh a bit, find all nodes in elements in element set UPPER and put 0.2 to X3 to test adjust
JuliaFEM.Preprocess.create_node_set_from_element_set!(mesh, :UPPER)
for nid in mesh.node_sets[:UPPER]
mesh.nodes[nid][3] += 0.2
end
upper = Problem(Elasticity, "UPPER", 3)
upper.elements = create_elements(mesh, "UPPER")
update!(upper, "youngs modulus", 288.0)
update!(upper, "poissons ratio", 1/3)
lower = Problem(Elasticity, "LOWER", 3)
lower.elements = create_elements(mesh, "LOWER")
update!(lower, "youngs modulus", 288.0)
update!(lower, "poissons ratio", 1/3)
bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement")
bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper, "displacement 3", 0.0)
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement")
bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower, "displacement 3", 0.0)
bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement")
bc_sym13.elements = [create_surface_elements(mesh, "LOWER_SYM13"); create_surface_elements(mesh, "UPPER_SYM13")]
update!(bc_sym13, "displacement 2", 0.0)
bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement")
bc_sym23.elements = [create_surface_elements(mesh, "LOWER_SYM23"); create_surface_elements(mesh, "UPPER_SYM23")]
update!(bc_sym23, "displacement 1", 0.0)
interface = Problem(Mortar, "LOWER_TO_UPPER", 3, "displacement")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
interface.elements = [interface_slave_elements; interface_master_elements]
append!(interface.assembly.removed_dofs, [1316, 1319, 1358, 1387, 1388, 1492, 1597, 1627])
interface.properties.linear_surface_elements = false
interface.properties.split_quadratic_slave_elements = false
interface.properties.split_quadratic_master_elements = false
interface.properties.adjust = true
interface.properties.dual_basis = true
solver = LinearSolver(upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface)
solver.xdmf = Xdmf("dl_lin_disp_results")
solver()
node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0)
node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0)
u3 = [u[3] for u in displacement]
maxabsu3 = maximum(abs(u3))
stdabsu3 = std(abs(u3))
info("tet10 block: max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-10)
end
@testset "patch test displacement + abaqus inp + tet10 + adjust" begin
mesh = abaqus_read_mesh(tet10_meshfile)
# modify mesh a bit, find all nodes in elements in element set UPPER and put 0.2 to X3 to test adjust
JuliaFEM.Preprocess.create_node_set_from_element_set!(mesh, :UPPER)
for nid in mesh.node_sets[:UPPER]
mesh.nodes[nid][3] += 0.2
end
upper = Problem(Elasticity, "UPPER", 3)
upper.elements = create_elements(mesh, "UPPER")
update!(upper, "youngs modulus", 288.0)
update!(upper, "poissons ratio", 1/3)
lower = Problem(Elasticity, "LOWER", 3)
lower.elements = create_elements(mesh, "LOWER")
update!(lower, "youngs modulus", 288.0)
update!(lower, "poissons ratio", 1/3)
bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement")
bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper, "displacement 3", 0.0)
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement")
bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower, "displacement 3", 0.0)
bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement")
bc_sym13.elements = [create_surface_elements(mesh, "LOWER_SYM13"); create_surface_elements(mesh, "UPPER_SYM13")]
update!(bc_sym13, "displacement 2", 0.0)
bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement")
bc_sym23.elements = [create_surface_elements(mesh, "LOWER_SYM23"); create_surface_elements(mesh, "UPPER_SYM23")]
update!(bc_sym23, "displacement 1", 0.0)
interface = Problem(Mortar, "LOWER_TO_UPPER", 3, "displacement")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
interface.elements = [interface_slave_elements; interface_master_elements]
removed_dofs = [1316, 1319, 1325, 1358, 1361, 1387, 1388, 1391, 1492, 1597, 1600, 1627, 1630, 1657]
append!(interface.assembly.removed_dofs, removed_dofs)
interface.properties.linear_surface_elements = false
interface.properties.split_quadratic_slave_elements = false
interface.properties.split_quadratic_master_elements = false
interface.properties.adjust = true
interface.properties.dual_basis = false
solver = LinearSolver(upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface)
solver.xdmf = Xdmf("sl_quad_disp_results")
solver()
node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0)
node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0)
u3 = [u[3] for u in displacement]
maxabsu3 = maximum(abs(u3))
stdabsu3 = std(abs(u3))
info("tet10 block: max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-6)
end
@testset "patch test displacement + abaqus inp + tet10 + adjust + dual basis + alpha=0.2" begin
mesh = abaqus_read_mesh(tet10_meshfile)
# modify mesh a bit, find all nodes in elements in element set UPPER and put 0.2 to X3 to test adjust
JuliaFEM.Preprocess.create_node_set_from_element_set!(mesh, :UPPER)
for nid in mesh.node_sets[:UPPER]
mesh.nodes[nid][3] += 0.2
end
upper = Problem(Elasticity, "UPPER", 3)
upper.elements = create_elements(mesh, "UPPER")
update!(upper, "youngs modulus", 288.0)
update!(upper, "poissons ratio", 1/3)
lower = Problem(Elasticity, "LOWER", 3)
lower.elements = create_elements(mesh, "LOWER")
update!(lower, "youngs modulus", 288.0)
update!(lower, "poissons ratio", 1/3)
bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement")
bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper, "displacement 3", 0.0)
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement")
bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower, "displacement 3", 0.0)
bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement")
bc_sym13.elements = [create_surface_elements(mesh, "LOWER_SYM13"); create_surface_elements(mesh, "UPPER_SYM13")]
update!(bc_sym13, "displacement 2", 0.0)
bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement")
bc_sym23.elements = [create_surface_elements(mesh, "LOWER_SYM23"); create_surface_elements(mesh, "UPPER_SYM23")]
update!(bc_sym23, "displacement 1", 0.0)
interface = Problem(Mortar, "LOWER_TO_UPPER", 3, "displacement")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
interface.elements = [interface_slave_elements; interface_master_elements]
removed_dofs = [1316, 1319, 1325, 1358, 1361, 1387, 1388, 1391, 1492, 1597, 1600, 1627, 1630, 1657]
append!(interface.assembly.removed_dofs, removed_dofs)
interface.properties.linear_surface_elements = false
interface.properties.split_quadratic_slave_elements = false
interface.properties.split_quadratic_master_elements = false
interface.properties.adjust = true
interface.properties.dual_basis = true
interface.properties.alpha = 0.2
solver = LinearSolver(upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface)
solver.xdmf = Xdmf("dl_quad_disp_results")
solver()
node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0)
node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0)
u3 = [u[3] for u in displacement]
maxabsu3 = maximum(abs(u3))
stdabsu3 = std(abs(u3))
info("tet10 block: max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-6)
end
test/testdata/test_problems_mortar_3d_tet10.inp → test/test_problems_mortar_3d/tet10.inp
Executable → Regular
View File
test/test_problems_mortar_3d/tet4.inp
+589
View File
@@ -0,0 +1,589 @@
**NSET COUNT = 116
*NODE
127, 1.00000, 0.00000, 0.75000
128, 1.00000, 0.00000, 1.00000
129, 0.75000, 0.00000, 1.00000
133, 0.75000, 0.00000, 0.75000
136, 1.00000, 0.00000, 0.50000
139, 0.75000, 0.00000, 0.50000
142, 0.50000, 0.00000, 1.00000
145, 0.50000, 0.00000, 0.75000
149, 0.50000, 0.00000, 0.50000
152, 0.25000, 0.00000, 1.00000
155, 0.25000, 0.00000, 0.75000
159, 0.25000, 0.00000, 0.50000
162, 0.00000, 0.00000, 1.00000
165, 0.00000, 0.00000, 0.75000
169, 0.00000, 0.00000, 0.50000
172, 0.75000, 1.00000, 1.00000
173, 1.00000, 1.00000, 1.00000
174, 1.00000, 1.00000, 0.75000
178, 0.75000, 1.00000, 0.75000
181, 1.00000, 1.00000, 0.50000
184, 0.75000, 1.00000, 0.50000
187, 0.50000, 1.00000, 1.00000
190, 0.50000, 1.00000, 0.75000
194, 0.50000, 1.00000, 0.50000
197, 0.25000, 1.00000, 1.00000
200, 0.25000, 1.00000, 0.75000
204, 0.25000, 1.00000, 0.50000
207, 0.00000, 1.00000, 1.00000
210, 0.00000, 1.00000, 0.75000
214, 0.00000, 1.00000, 0.50000
217, 1.00000, 0.75000, 1.00000
220, 1.00000, 0.75000, 0.75000
224, 1.00000, 0.75000, 0.50000
227, 1.00000, 0.50000, 1.00000
230, 1.00000, 0.50000, 0.75000
234, 1.00000, 0.50000, 0.50000
237, 1.00000, 0.25000, 1.00000
240, 1.00000, 0.25000, 0.75000
244, 1.00000, 0.25000, 0.50000
252, 0.00000, 0.75000, 1.00000
255, 0.00000, 0.75000, 0.75000
259, 0.00000, 0.75000, 0.50000
262, 0.00000, 0.50000, 1.00000
265, 0.00000, 0.50000, 0.75000
269, 0.00000, 0.50000, 0.50000
272, 0.00000, 0.25000, 1.00000
275, 0.00000, 0.25000, 0.75000
279, 0.00000, 0.25000, 0.50000
288, 0.75000, 0.75000, 1.00000
292, 0.50000, 0.75000, 1.00000
296, 0.25000, 0.75000, 1.00000
302, 0.75000, 0.50000, 1.00000
306, 0.50000, 0.50000, 1.00000
310, 0.25000, 0.50000, 1.00000
316, 0.75000, 0.25000, 1.00000
320, 0.50000, 0.25000, 1.00000
324, 0.25000, 0.25000, 1.00000
337, 0.75000, 0.75000, 0.50000
341, 0.50000, 0.75000, 0.50000
345, 0.25000, 0.75000, 0.50000
351, 0.75000, 0.50000, 0.50000
355, 0.50000, 0.50000, 0.50000
359, 0.25000, 0.50000, 0.50000
365, 0.75000, 0.25000, 0.50000
369, 0.50000, 0.25000, 0.50000
373, 0.25000, 0.25000, 0.50000
438, 1.00000, 0.00000, 0.25000
439, 1.00000, 0.00000, 0.50000
440, 0.66667, 0.00000, 0.50000
444, 0.66667, 0.00000, 0.25000
447, 1.00000, 0.00000, 0.00000
450, 0.66667, 0.00000, 0.00000
453, 0.33333, 0.00000, 0.50000
456, 0.33333, 0.00000, 0.25000
460, 0.33333, 0.00000, 0.00000
463, 0.00000, 0.00000, 0.50000
466, 0.00000, 0.00000, 0.25000
470, 0.00000, 0.00000, 0.00000
473, 0.66667, 1.00000, 0.50000
474, 1.00000, 1.00000, 0.50000
475, 1.00000, 1.00000, 0.25000
479, 0.66667, 1.00000, 0.25000
482, 1.00000, 1.00000, 0.00000
485, 0.66667, 1.00000, 0.00000
488, 0.33333, 1.00000, 0.50000
491, 0.33333, 1.00000, 0.25000
495, 0.33333, 1.00000, 0.00000
498, 0.00000, 1.00000, 0.50000
501, 0.00000, 1.00000, 0.25000
505, 0.00000, 1.00000, 0.00000
508, 1.00000, 0.66667, 0.50000
511, 1.00000, 0.66667, 0.25000
515, 1.00000, 0.66667, 0.00000
518, 1.00000, 0.33333, 0.50000
521, 1.00000, 0.33333, 0.25000
525, 1.00000, 0.33333, 0.00000
533, 0.00000, 0.66667, 0.50000
536, 0.00000, 0.66667, 0.25000
540, 0.00000, 0.66667, 0.00000
543, 0.00000, 0.33333, 0.50000
546, 0.00000, 0.33333, 0.25000
550, 0.00000, 0.33333, 0.00000
559, 0.66667, 0.66667, 0.50000
563, 0.33333, 0.66667, 0.50000
569, 0.66667, 0.33333, 0.50000
573, 0.33333, 0.33333, 0.50000
584, 0.66667, 0.66667, 0.00000
588, 0.33333, 0.66667, 0.00000
594, 0.66667, 0.33333, 0.00000
598, 0.33333, 0.33333, 0.00000
608, 0.66006, 0.47128, 0.70078
609, 0.19762, 0.64074, 0.81889
610, 0.65771, 0.82813, 0.74829
611, 0.59993, 0.15953, 0.76106
612, 0.40915, 0.16311, 0.74937
613, 0.50000, 0.50000, 0.25000
**
**ELSET COUNT = 172
**HWCOLOR COMP 54 0
*ELEMENT, TYPE=C3D4, ELSET=UPPER
570, 252, 262, 255, 609
571, 252, 310, 262, 609
572, 252, 296, 310, 609
573, 200, 252, 255, 609
574, 296, 306, 310, 609
575, 262, 310, 265, 609
576, 259, 265, 269, 359
577, 184, 337, 220, 610
578, 259, 345, 265, 359
579, 288, 608, 292, 610
580, 187, 288, 292, 610
581, 172, 288, 187, 610
582, 178, 172, 187, 610
583, 178, 220, 172, 610
584, 172, 220, 288, 610
585, 230, 240, 237, 316
586, 178, 184, 220, 610
587, 178, 190, 184, 610
588, 178, 187, 190, 610
589, 187, 292, 190, 610
590, 190, 292, 609, 610
591, 190, 609, 341, 610
592, 190, 341, 194, 610
593, 184, 190, 194, 610
594, 184, 194, 337, 610
595, 194, 341, 337, 610
596, 337, 341, 608, 610
597, 230, 337, 608, 610
598, 230, 608, 288, 610
599, 220, 230, 288, 610
600, 220, 337, 230, 610
601, 341, 609, 608, 610
602, 292, 608, 609, 610
603, 306, 316, 320, 611
604, 306, 608, 316, 611
605, 365, 608, 369, 611
606, 240, 608, 365, 611
607, 240, 316, 608, 611
608, 240, 133, 316, 611
609, 240, 365, 133, 611
610, 365, 139, 133, 611
611, 365, 369, 139, 611
612, 369, 149, 139, 611
613, 149, 145, 611, 612
614, 240, 127, 128, 316
615, 320, 129, 142, 611
616, 316, 129, 320, 611
617, 316, 133, 129, 611
618, 129, 133, 142, 611
619, 142, 133, 145, 611
620, 139, 145, 133, 611
621, 139, 149, 145, 611
622, 310, 320, 324, 612
623, 306, 320, 310, 612
624, 306, 611, 320, 612
625, 306, 608, 611, 612
626, 306, 609, 608, 612
627, 306, 310, 609, 612
628, 265, 609, 310, 612
629, 265, 310, 324, 612
630, 265, 324, 275, 612
631, 265, 275, 373, 612
632, 265, 373, 359, 612
633, 265, 359, 609, 612
634, 355, 609, 359, 612
635, 355, 608, 609, 612
636, 355, 369, 608, 612
637, 369, 611, 608, 612
638, 214, 255, 259, 345
639, 142, 611, 145, 612
640, 320, 611, 142, 612
641, 320, 142, 324, 612
642, 324, 142, 152, 612
643, 324, 152, 155, 612
644, 275, 324, 155, 612
645, 275, 155, 373, 612
646, 373, 155, 159, 612
647, 373, 159, 149, 612
648, 369, 373, 149, 612
649, 359, 373, 369, 612
650, 355, 359, 369, 612
651, 369, 149, 611, 612
652, 149, 155, 145, 612
653, 149, 159, 155, 612
654, 152, 145, 155, 612
655, 142, 145, 152, 612
656, 237, 240, 128, 316
657, 244, 136, 127, 365
658, 240, 365, 127, 133
659, 240, 351, 244, 365
660, 220, 224, 230, 337
661, 227, 288, 230, 302
662, 200, 207, 252, 296
663, 240, 127, 316, 133
664, 230, 337, 234, 351
665, 275, 162, 165, 155
666, 172, 220, 217, 288
667, 200, 252, 210, 255
668, 190, 292, 197, 296
669, 230, 234, 240, 351
670, 204, 210, 214, 345
671, 210, 255, 214, 345
672, 275, 165, 373, 155
673, 172, 174, 217, 220
674, 200, 204, 341, 345
675, 240, 244, 127, 365
676, 220, 230, 227, 288
677, 272, 275, 324, 155
678, 272, 152, 162, 155
679, 200, 255, 210, 345
680, 190, 200, 194, 341
681, 197, 207, 200, 296
682, 190, 197, 200, 296
683, 279, 165, 169, 159
684, 200, 207, 210, 252
685, 224, 234, 230, 337
686, 187, 197, 190, 292
687, 275, 165, 279, 373
688, 265, 275, 269, 373
689, 272, 324, 152, 155
690, 262, 265, 310, 324
691, 265, 272, 275, 324
692, 184, 224, 220, 337
693, 272, 162, 275, 155
694, 262, 272, 265, 324
695, 178, 174, 172, 220
696, 136, 365, 139, 133
697, 227, 230, 237, 302
698, 178, 184, 181, 220
699, 172, 174, 173, 217
700, 194, 200, 204, 341
701, 234, 244, 240, 351
702, 200, 210, 204, 345
703, 230, 237, 302, 316
704, 128, 316, 127, 129
705, 217, 220, 227, 288
706, 127, 365, 136, 133
707, 265, 269, 359, 373
708, 279, 373, 165, 159
709, 165, 159, 373, 155
710, 269, 275, 279, 373
711, 127, 129, 316, 133
712, 181, 220, 184, 224
713, 178, 181, 174, 220
714, 292, 302, 306, 608
715, 288, 302, 292, 608
716, 230, 302, 288, 608
717, 230, 316, 302, 608
718, 230, 240, 316, 608
719, 230, 351, 240, 608
720, 230, 337, 351, 608
721, 302, 316, 306, 608
722, 240, 351, 365, 608
723, 351, 355, 365, 608
724, 341, 355, 351, 608
725, 337, 341, 351, 608
726, 355, 369, 365, 608
727, 341, 355, 608, 609
728, 292, 608, 306, 609
729, 292, 306, 296, 609
730, 190, 292, 296, 609
731, 190, 200, 341, 609
732, 190, 296, 200, 609
733, 200, 296, 252, 609
734, 200, 255, 345, 609
735, 200, 345, 341, 609
736, 341, 345, 355, 609
737, 345, 359, 355, 609
738, 265, 359, 345, 609
739, 255, 265, 345, 609
740, 255, 265, 259, 345
741, 255, 262, 265, 609
**
**ELSET COUNT = 92
**HWCOLOR COMP 1 0
*ELEMENT, TYPE=C3D4, ELSET=LOWER
742, 456, 460, 466, 598
743, 444, 450, 456, 598
744, 444, 525, 447, 594
745, 536, 540, 588, 598
746, 453, 456, 463, 543
747, 533, 543, 536, 573
748, 444, 456, 453, 573
749, 466, 550, 546, 598
750, 536, 543, 546, 573
751, 491, 536, 495, 588
752, 438, 447, 444, 525
753, 444, 594, 450, 598
754, 444, 521, 525, 594
755, 438, 444, 440, 521
756, 540, 546, 550, 598
757, 453, 456, 543, 573
758, 488, 491, 559, 563
759, 505, 536, 540, 588
760, 440, 444, 453, 573
761, 491, 495, 584, 588
762, 444, 447, 450, 594
763, 485, 491, 495, 584
764, 440, 569, 444, 573
765, 495, 501, 505, 536
766, 473, 479, 511, 559
767, 438, 440, 439, 521
768, 491, 533, 536, 563
769, 444, 569, 521, 594
770, 473, 479, 475, 511
771, 491, 498, 501, 533
772, 495, 536, 505, 588
773, 438, 444, 521, 525
774, 473, 511, 508, 559
775, 440, 521, 444, 569
776, 511, 515, 521, 584
777, 491, 498, 533, 563
778, 515, 525, 521, 584
779, 460, 550, 466, 598
780, 439, 521, 440, 569
781, 456, 463, 543, 546
782, 479, 482, 511, 515
783, 533, 536, 563, 573
784, 479, 515, 511, 584
785, 450, 460, 456, 598
786, 521, 569, 559, 594
787, 559, 594, 569, 613
788, 559, 584, 594, 613
789, 521, 559, 584, 594
790, 536, 546, 540, 598
791, 511, 521, 559, 584
792, 479, 559, 491, 584
793, 491, 559, 563, 584
794, 511, 521, 518, 559
795, 491, 584, 563, 588
796, 563, 588, 584, 613
797, 559, 563, 584, 613
798, 559, 569, 563, 613
799, 563, 569, 573, 613
800, 536, 563, 573, 588
801, 563, 573, 588, 613
802, 573, 598, 588, 613
803, 491, 501, 495, 536
804, 536, 573, 546, 598
805, 536, 588, 573, 598
806, 588, 598, 594, 613
807, 584, 588, 594, 613
808, 444, 569, 594, 598
809, 521, 584, 525, 594
810, 444, 573, 569, 598
811, 569, 598, 573, 613
812, 569, 594, 598, 613
813, 444, 456, 573, 598
814, 439, 518, 521, 569
815, 456, 546, 543, 573
816, 479, 485, 482, 515
817, 508, 511, 518, 559
818, 456, 546, 573, 598
819, 488, 498, 491, 563
820, 479, 485, 515, 584
821, 491, 563, 536, 588
822, 473, 488, 479, 559
823, 479, 491, 485, 584
824, 491, 533, 501, 536
825, 479, 511, 559, 584
826, 479, 488, 491, 559
827, 475, 479, 482, 511
828, 518, 559, 521, 569
829, 473, 475, 508, 511
830, 473, 475, 474, 508
831, 460, 470, 466, 550
832, 456, 466, 463, 546
833, 456, 466, 546, 598
**
**ELSET COUNT = 18
*SURFACE,TYPE=ELEMENT,NAME=LOWER_TO_UPPER
799,S1
798,S1
764,S2
814,S2
747,S2
828,S2
819,S2
822,S2
777,S3
746,S4
757,S4
760,S4
780,S4
783,S4
817,S4
758,S4
774,S4
830,S4
**
**ELSET COUNT = 18
*SURFACE,TYPE=ELEMENT,NAME=LOWER_BOTTOM
806,S1
807,S1
831,S2
779,S2
785,S2
778,S2
753,S3
762,S3
744,S3
745,S3
809,S3
761,S3
820,S3
816,S3
756,S4
759,S4
772,S4
763,S4
**
**ELSET COUNT = 12
*SURFACE,TYPE=ELEMENT,NAME=LOWER_SYM13
831,S1
742,S1
832,S1
746,S1
785,S1
743,S1
748,S1
760,S1
762,S1
752,S1
755,S1
767,S1
**
**ELSET COUNT = 12
*SURFACE,TYPE=ELEMENT,NAME=LOWER_SYM23
749,S1
756,S1
790,S1
750,S1
747,S1
759,S1
831,S3
832,S3
781,S3
765,S3
824,S3
771,S3
**
**ELSET COUNT = 32
*SURFACE,TYPE=ELEMENT,NAME=UPPER_TOP
678,S1
689,S1
642,S1
641,S1
615,S1
616,S1
622,S1
623,S1
603,S1
721,S1
571,S1
572,S1
574,S1
729,S1
714,S1
715,S1
580,S1
581,S1
704,S2
694,S2
661,S2
681,S2
686,S2
703,S3
662,S3
668,S3
656,S4
690,S4
697,S4
705,S4
666,S4
699,S4
**
**ELSET COUNT = 32
*SURFACE,TYPE=ELEMENT,NAME=UPPER_TO_LOWER
647,S1
648,S1
612,S1
611,S1
696,S1
649,S1
650,S1
726,S1
723,S1
737,S1
736,S1
724,S1
725,S1
595,S1
594,S1
708,S2
657,S2
701,S2
578,S2
685,S2
692,S2
707,S3
659,S3
664,S3
674,S3
683,S4
710,S4
576,S4
638,S4
670,S4
700,S4
712,S4
**
**ELSET COUNT = 16
*SURFACE,TYPE=ELEMENT,NAME=UPPER_SYM13
653,S1
652,S1
654,S1
655,S1
621,S1
620,S1
619,S1
618,S1
709,S2
711,S2
683,S3
665,S3
678,S3
696,S4
706,S4
704,S4
**
**ELSET COUNT = 16
*SURFACE,TYPE=ELEMENT,NAME=UPPER_SYM23
683,S1
687,S1
665,S1
693,S1
710,S1
688,S1
691,S1
694,S1
576,S1
740,S1
741,S1
570,S1
638,S1
671,S1
667,S3
684,S3
**
**Property Definitions
**
*SOLID SECTION, ELSET=UPPER, MATERIAL=Def_Material
*SOLID SECTION, ELSET=LOWER, MATERIAL=Def_Material
**
**Material Definitions
**
**Material:Def_Material
*MATERIAL,NAME=Def_Material
*ELASTIC,TYPE=ISO
2.08000e+005,3.00000e-001
*DENSITY
7.80000e-009,
*SPECIFIC HEAT
5.00000e-001
*CONDUCTIVITY
4.98100e-002
**
test/test_problems_mortar_3d_lowlevel.jl
+196
View File
@@ -0,0 +1,196 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Postprocess
using JuliaFEM.Testing
using JuliaFEM.Abaqus: create_surface_elements
@testset "forget to add elements to problem" begin
X = Dict(
1 => [0.0, 0.0, 0.0],
2 => [1.0, 0.0, 0.0],
3 => [0.0, 1.0, 0.0],
4 => [-0.25, 0.50, 0.00],
5 => [0.50, -0.25, 0.00],
6 => [0.75, 0.75, 0.00])
s = Element(Tri3, [1, 2, 3])
m = Element(Tri3, [4, 5, 6])
update!([s, m], "geometry", X)
update!(s, "master elements", [m])
p = Problem(Mortar, "two elements", 1, "temperature")
initialize!(p)
assemble!(p)
@test true
end
""" Calculate mortar projection matrix P = D^-1*M from mortar assembly. """
function calculate_mortar_projection_matrix(problem::Problem{Mortar}, ndim::Int)
C1 = sparse(problem.assembly.C1, ndim, ndim)
C2 = sparse(problem.assembly.C2, ndim, ndim)
@assert nnz(sparse(problem.assembly.K)) == 0
@assert nnz(sparse(problem.assembly.D)) == 0
@assert nnz(sparse(problem.assembly.Kg)) == 0
@assert nnz(sparse(problem.assembly.fg)) == 0
@assert nnz(sparse(problem.assembly.f)) == 0
@assert nnz(sparse(problem.assembly.g)) == 0
@assert C1 == C2
#@assert problem.properties.dual_basis == true
@assert problem.properties.adjust == false
S = get_nonzero_rows(C2)
M = setdiff(get_nonzero_columns(C2), S)
# Construct matrix P = D^-1*M
D_ = C2[S,S]
M_ = -C2[S,M]
#=
P = nothing
if !isdiag(D_)
warn("D is not diagonal, is dual basis used? This might take a long time.")
P = ldltfact(1/2*(D_ + D_')) \ M_
else
P = D_ \ M_
end
=#
P = lufact(D_) \ full(M_)
return S, M, P
end
@testset "two linear element clipping, calculation of projection matrix P for standard and dual basis" begin
X = Dict(
1 => [0.0, 0.0, 0.0],
2 => [1.0, 0.0, 0.0],
3 => [0.0, 1.0, 0.0],
4 => [-0.25, 0.50, 0.00],
5 => [0.50, -0.25, 0.00],
6 => [0.75, 0.75, 0.00])
s = Element(Tri3, [1, 2, 3])
m = Element(Tri3, [4, 5, 6])
update!([s, m], "geometry", X)
update!(s, "master elements", [m])
p = Problem(Mortar, "two elements", 1, "temperature")
p.properties.dual_basis = false
p.elements = [s; m]
initialize!(p)
assemble!(p)
C1 = sparse(p.assembly.C1)
C2 = sparse(p.assembly.C2)
D = sparse(p.assembly.D)
@test length(D) == 0
@test C1 == C2
S, M, P = calculate_mortar_projection_matrix(p, 6)
@test S == [1, 2, 3]
@test M == [4, 5, 6]
# visually inspected to be ok result
P_expected = 1/15*[9 9 -3; -7 13 9; 13 -7 9]
@test isapprox(P, P_expected)
um = [7.5, 15.0, 22.5]
@test isapprox(P*um, [9.0, 23.0, 13.0])
empty!(p.assembly)
p.properties.dual_basis = true
assemble!(p)
C1 = sparse(p.assembly.C1)
C2 = sparse(p.assembly.C2)
D = sparse(p.assembly.D)
@test length(D) == 0
@test C1 == C2
S, M, P = calculate_mortar_projection_matrix(p, 6)
@test S == [1, 2, 3]
@test M == [4, 5, 6]
@test isapprox(P, P_expected)
um = [7.5, 15.0, 22.5]
@test isapprox(P*um, [9.0, 23.0, 13.0])
end
@testset "two quadratic element clipping, calculation of projection matrix P for standard basis" begin
X = Dict(
1 => [0.0, 0.0, 0.0],
2 => [1.0, 0.0, 0.0],
3 => [0.0, 1.0, 0.0],
7 => [-0.25, 0.50, 0.00],
8 => [0.50, -0.25, 0.00],
9 => [0.75, 0.75, 0.00])
# middle nodes
X[4] = 1/2*(X[1] + X[2])
X[5] = 1/2*(X[2] + X[3])
X[6] = 1/2*(X[3] + X[1])
X[10] = 1/2*(X[7] + X[8])
X[11] = 1/2*(X[8] + X[9])
X[12] = 1/2*(X[9] + X[7])
s = Element(Tri6, [1, 2, 3, 4, 5, 6])
m = Element(Tri6, [7, 8, 9, 10, 11, 12])
update!([s, m], "geometry", X)
update!(s, "master elements", [m])
p = Problem(Mortar, "two elements", 1, "temperature")
p.properties.dual_basis = false
p.properties.alpha = 0.2
p.elements = [s; m]
initialize!(p)
assemble!(p)
C1 = sparse(p.assembly.C1)
C2 = sparse(p.assembly.C2)
D = sparse(p.assembly.D)
@test length(D) == 0
@test C1 == C2
S, M, P = calculate_mortar_projection_matrix(p, 12)
@test S == [1, 2, 3, 4, 5, 6]
@test M == [7, 8, 9, 10, 11, 12]
println(full(P))
# visually inspected to be ok result
P_expected = 1/675*[81 81 189 972 -324 -324; 609 429 81 -1092 1404 -756; 429 609 81 -1092 -756 1404; -39 231 -81 132 396 36; -81 -81 81 108 324 324; 231 -39 -81 132 36 396]
@test isapprox(P, P_expected)
um = 15/2*[1, 2, 3]
um = [um[1], um[2], um[3], 0.5*(um[1]+um[2]), 0.5*(um[2]+um[3]), 0.5*(um[3]+um[1])]
us = P*um
@test isapprox(us, [9.0, 23.0, 13.0, 16.0, 18.0, 11.0])
end
@testset "two quadratic element clipping, calculation of projection matrix P for dual lagrange basis" begin
X = Dict(
1 => [0.0, 0.0, 0.0],
2 => [1.0, 0.0, 0.0],
3 => [0.0, 1.0, 0.0],
7 => [-0.25, 0.50, 0.00],
8 => [0.50, -0.25, 0.00],
9 => [0.75, 0.75, 0.00])
# middle nodes
X[4] = 1/2*(X[1] + X[2])
X[5] = 1/2*(X[2] + X[3])
X[6] = 1/2*(X[3] + X[1])
X[10] = 1/2*(X[7] + X[8])
X[11] = 1/2*(X[8] + X[9])
X[12] = 1/2*(X[9] + X[7])
s = Element(Tri6, [1, 2, 3, 4, 5, 6])
m = Element(Tri6, [7, 8, 9, 10, 11, 12])
update!([s, m], "geometry", X)
update!(s, "master elements", [m])
p = Problem(Mortar, "two elements", 1, "temperature")
p.properties.dual_basis = true
p.properties.alpha = 0.2
p.elements = [s; m]
initialize!(p)
assemble!(p)
C1 = sparse(p.assembly.C1)
C2 = sparse(p.assembly.C2)
D = sparse(p.assembly.D)
@test length(D) == 0
@test C1 == C2
S, M, P = calculate_mortar_projection_matrix(p, 12)
P_expected = 1/675*[81 81 189 972 -324 -324; 609 429 81 -1092 1404 -756; 429 609 81 -1092 -756 1404; -39 231 -81 132 396 36; -81 -81 81 108 324 324; 231 -39 -81 132 36 396]
@test isapprox(P, P_expected)
um = 15/2*[1, 2, 3]
um = [um[1], um[2], um[3], 0.5*(um[1]+um[2]), 0.5*(um[2]+um[3]), 0.5*(um[3]+um[1])]
us = P*um
@test isapprox(us, [9.0, 23.0, 13.0, 16.0, 18.0, 11.0])
end