mortar tests pass now also

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")

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

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
 

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
 
-

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