diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 6ab0373..03a8f9c 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -84,7 +84,7 @@ include("solvers.jl") include("directsolver.jl") # parallel sparse direct solver for non-linear problems ### MORTAR STUFF ### -#include("mortar.jl") # mortar projection +include("mortar.jl") # mortar projection # PRE AND POSTPROCESS include("xdmf.jl") diff --git a/src/dirichlet.jl b/src/dirichlet.jl index 1c2312c..e7749a7 100644 --- a/src/dirichlet.jl +++ b/src/dirichlet.jl @@ -7,7 +7,7 @@ function DirichletProblem(parent_field_name, parent_field_dim, dim=1, elements=[ return BoundaryProblem{DirichletProblem}(parent_field_name, parent_field_dim, dim, elements) end -function assemble!{E}(assembly::Assembly, problem::BoundaryProblem{DirichletProblem}, element::Element{E}, time::Number) +function assemble!(assembly::Assembly, problem::BoundaryProblem{DirichletProblem}, element::Element, time::Number) # get dimension and name of PARENT field field_dim = problem.parent_field_dim diff --git a/src/elements.jl b/src/elements.jl index fdb1028..00c7b72 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -137,7 +137,11 @@ function get_basis{E}(element::Element{E}, ip::IntegrationPoint) return get_basis(E, ip.xi) end -function call{E}(element::Element{E}, xi::VecOrIP, time::Float64=0) +function get_basis{E}(element::Element{E}, xi::Vector{Float64}) + return get_basis(E, xi) +end + +function call{E}(element::Element{E}, xi::VecOrIP, time::Float64=0.0) return get_basis(element, xi) end diff --git a/src/mortar.jl b/src/mortar.jl index 15cd9e4..30a5a8d 100644 --- a/src/mortar.jl +++ b/src/mortar.jl @@ -6,7 +6,7 @@ """ Find projection from slave nodes to master element, i.e. find xi2 from master element corresponding to the xi1. """ -function project_from_slave_to_master(slave::Element, master::Element, xi1::Vector, time::Float64=0.0; max_iterations=5, tol=1.0e-9) +function project_from_slave_to_master{S,M}(slave::Element{S}, master::Element{M}, xi1::Vector, time::Float64=0.0; max_iterations=5, tol=1.0e-9) # slave_basis = get_basis(slave) # slave side geometry and normal direction at xi1 @@ -14,17 +14,19 @@ function project_from_slave_to_master(slave::Element, master::Element, xi1::Vect N1 = slave("nodal ntsys", xi1, time)[:,1] # master side geometry at xi2 - master_basis = master.basis.data.basis - master_dbasis = master.basis.data.dbasis + #master_basis = master.basis.data.basis + #master_dbasis = master.basis.data.dbasis + master_basis(xi) = get_basis(M, [xi]) + master_dbasis(xi) = get_dbasis(M, [xi]) master_geometry = master("geometry")(time) function X2(xi2) - N = master_basis([xi2]) + N = master_basis(xi2) return sum([N[i]*master_geometry[i] for i=1:length(N)]) end function dX2(xi2) - dN = master_dbasis([xi2]) + dN = master_dbasis(xi2) return sum([dN[i]*master_geometry[i] for i=1:length(dN)]) end @@ -51,33 +53,35 @@ end """ Find projection from master surface to slave point, i.e. find xi1 from slave element corresponding to the xi2. """ -function project_from_master_to_slave(slave::Element, master::Element, xi2::Vector, time::Float64=0.0; max_iterations=5, tol=1.0e-9) +function project_from_master_to_slave{S,M}(slave::Element{S}, master::Element{M}, xi2::Vector, time::Float64=0.0; max_iterations=5, tol=1.0e-9) # slave_basis = get_basis(slave) # slave side geometry and normal direction at xi1 slave_geometry = slave("geometry")(time) slave_normals = slave("nodal ntsys")(time) - slave_basis = slave.basis.data.basis - slave_dbasis = slave.basis.data.dbasis + #slave_basis = slave.basis.data.basis + #slave_dbasis = slave.basis.data.dbasis + slave_basis(xi) = get_basis(S, [xi]) + slave_dbasis(xi) = get_dbasis(S, [xi]) function X1(xi1) - N = slave_basis([xi1]) + N = slave_basis(xi1) return sum([N[i]*slave_geometry[i] for i=1:length(N)]) end function dX1(xi1) - dN = slave_dbasis([xi1]) + dN = slave_dbasis(xi1) return sum([dN[i]*slave_geometry[i] for i=1:length(dN)]) end function N1(xi1) - N = slave_basis([xi1]) + N = slave_basis(xi1) return sum([N[i]*slave_normals[i] for i=1:length(N)])[:,1] end function dN1(xi1) - dN = slave_dbasis([xi1]) + dN = slave_dbasis(xi1) return sum([dN[i]*slave_normals[i] for i=1:length(dN)])[:,1] end @@ -112,6 +116,7 @@ function project_from_master_to_slave(slave::Element, master::Element, xi2::Vect for i=1:max_iterations dxi1 = -R(xi1) / dR(xi1) xi1 += dxi1 + #info("dxi1 = $dxi1, xi1 = $xi1, norm(dxi1) = $(norm(dxi1))") if norm(dxi1) < tol return Float64[xi1] end @@ -119,30 +124,6 @@ function project_from_master_to_slave(slave::Element, master::Element, xi2::Vect error("find projection from master to slave: did not converge") end - -### Mortar equations - -abstract MortarEquation <: Equation - -""" Mortar boundary condition element for 2-dimensional problem, 2 node line segment. """ -type MBC2D2 <: MortarEquation - element :: Seg2 - integration_points :: Vector{IntegrationPoint} -end - -function Base.size(equation::MBC2D2) - return (1, 2) -end - -function Base.convert(::Type{MortarEquation}, element::Seg2) - integration_points = get_integration_points(element, Val{3}) - if !haskey(element, "reaction force") - element["reaction force"] = (0.0 => Vector{Float64}[]) - end - MBC2D2(element, integration_points) -end - - ### Mortar problem """ @@ -152,30 +133,23 @@ node_csys coordinate system in node, normal + tangent + "binormal" in 3d 3x3 matrix, in 2d 2x2 matrix, respectively """ -type MortarProblem <: BoundaryProblem - unknown_field_name :: ASCIIString - unknown_field_dimension :: Int - equations :: Vector{MortarEquation} -end +abstract MortarProblem <: AbstractProblem -function MortarProblem(unknown_field_name, unknown_field_dimension::Int=1) - MortarProblem(unknown_field_name, unknown_field_dimension, []) +function MortarProblem(parent_field_name, parent_field_dim, dim=1, elements=[]) + return BoundaryProblem{MortarProblem}(parent_field_name, parent_field_dim, dim, elements) end # Mortar assembly -function assemble!(assembly::Assembly, equation::MortarEquation, time::Number=0.0, problem=nothing) - isa(problem, Void) && error("Mortar boundary problem needs problem to be defined") - field_dim = problem.unknown_field_dimension - field_name = problem.unknown_field_name +function assemble!(assembly::Assembly, problem::BoundaryProblem{MortarProblem}, slave_element::Element, time::Number) + + # get dimension and name of PARENT field + field_dim = problem.parent_field_dim + field_name = problem.parent_field_name - slave_element = get_element(equation) slave_dofs = get_gdofs(slave_element, field_dim) - slave_basis = get_basis(slave_element) - detJ = det(slave_basis) for master_element in slave_element["master elements"] - master_dofs = get_gdofs(master_element, field_dim) xi1a = project_from_master_to_slave(slave_element, master_element, [-1.0]) xi1b = project_from_master_to_slave(slave_element, master_element, [ 1.0]) xi1 = clamp([xi1a xi1b], -1.0, 1.0) @@ -184,9 +158,9 @@ function assemble!(assembly::Assembly, equation::MortarEquation, time::Number=0. warn("No contribution") continue # no contribution end - master_basis = get_basis(master_element) - for ip in get_integration_points(equation) - w = ip.weight*detJ(ip)*l + master_dofs = get_gdofs(master_element, field_dim) + for ip in get_integration_points(slave_element) + w = ip.weight*det(slave_element, ip, time)*l # integration point on slave side segment xi_gauss = 1/2*(1-ip.xi)*xi1[1] + 1/2*(1+ip.xi)*xi1[2] @@ -194,8 +168,8 @@ function assemble!(assembly::Assembly, equation::MortarEquation, time::Number=0. xi_projected = project_from_slave_to_master(slave_element, master_element, xi_gauss) # add contribution to left hand side - N1 = slave_basis(xi_gauss, time) - N2 = master_basis(xi_projected, time) + N1 = slave_element(xi_gauss, time) + N2 = master_element(xi_projected, time) S = w*N1'*N1 M = w*N1'*N2 for i=1:field_dim @@ -209,4 +183,3 @@ function assemble!(assembly::Assembly, equation::MortarEquation, time::Number=0. end end - diff --git a/test/test_mortar.jl b/test/test_mortar.jl index 82a1cf1..0075fd0 100644 --- a/test/test_mortar.jl +++ b/test/test_mortar.jl @@ -6,9 +6,9 @@ module MortarTests using JuliaFEM using JuliaFEM.Test -using JuliaFEM: Seg2, MortarProblem, MortarEquation, MortarElement, Assembly, assemble!, Element -using JuliaFEM: get_basis, grad, project_from_slave_to_master, project_from_master_to_slave, Quad4 +using JuliaFEM: Element, Seg2, Quad4, MortarProblem, Assembly, assemble! using JuliaFEM: PlaneStressElasticityProblem, DirichletProblem, DirectSolver +using JuliaFEM: project_from_slave_to_master, project_from_master_to_slave function get_test_2d_model() # this is hand calculated and given as an example in my thesis @@ -40,7 +40,7 @@ function get_test_2d_model() end -function test_calc_flat_2d_projection() +function test_calc_flat_2d_projection_slave_to_master() slaves, masters = get_test_2d_model() slave1, slave2 = slaves master1, master2 = masters @@ -50,15 +50,21 @@ function test_calc_flat_2d_projection() xi2b = project_from_slave_to_master(slave1, master1, [1.0]) @test xi2b == [ 0.2] - X2 = get_basis(master1)("geometry", xi2b) + X2 = master1("geometry", xi2b, 0.0) @test X2 == [3/4, 1.0] +end +function test_calc_flat_2d_projection_master_to_slave() + slaves, masters = get_test_2d_model() + slave1, slave2 = slaves + master1, master2 = masters xi1a = project_from_master_to_slave(slave1, master1, [-1.0]) @test xi1a == [-1.0] xi1b = project_from_master_to_slave(slave1, master1, [1.0]) - X1 = get_basis(slave1)("geometry", xi1b) + X1 = slave1("geometry", xi1b, 0.0) @test X1 == [5/4, 1.0] end +#test_calc_flat_2d_projection_master_to_slave() function test_calc_flat_2d_projection_rotated() master1 = Seg2([3, 4]) @@ -80,7 +86,6 @@ function test_calc_flat_2d_projection_rotated() info("xi = $xi") @test xi == [-1.0] - end function test_create_flat_2d_assembly() @@ -102,7 +107,7 @@ function test_create_flat_2d_assembly() info("creating assembly") assembly = Assembly() - assemble!(assembly, problem.equations[1], 0.0, problem) + assemble!(assembly, problem, slave1, 0.0) B = round(full(assembly.stiffness_matrix, 12, 12), 6) info("size of B = $(size(B))") info("B matrix in first slave element = \n$(B[10:11,:])") @@ -120,7 +125,7 @@ function test_create_flat_2d_assembly() M3 = [8, 9] B_expected[S3,S3] += [9/100 27/200; 27/200 39/100] B_expected[S3,M3] -= [3/20 3/40; 9/40 3/10] - assemble!(assembly, problem.equations[2], 0.0, problem) + assemble!(assembly, problem, slave2, 0.0) B = full(assembly.stiffness_matrix) info("size of B = $(size(B))") info("B matrix in second slave element = \n$(B[11:12,:])") @@ -128,6 +133,7 @@ function test_create_flat_2d_assembly() @test isapprox(B, B_expected) end +#test_create_flat_2d_assembly() function test_2d_mortar_multiple_bodies_multiple_dirichlet_bc() N = Vector[ @@ -205,6 +211,7 @@ function test_2d_mortar_multiple_bodies_multiple_dirichlet_bc() # code aster verification, two_elements.comm @test isapprox(disp, [3.17431158889468E-02, -2.77183037855653E-01]) end +test_2d_mortar_multiple_bodies_multiple_dirichlet_bc() function test_2d_mortar_three_bodies_shared_nodes() @@ -325,6 +332,4 @@ function test_2d_mortar_three_bodies_shared_nodes() end - - end