Merge branch 'master' of git://github.com/JuliaFEM/JuliaFEM.jl into HEAD

src/assembly.jl

@@ -61,7 +61,9 @@ function assemble(problem::AllProblems, time::Real, nchunks=10)
     for (j, elrange) in enumerate(slices)
         sub_assembly = assemble(problem, elrange, time)
         append!(assembly, sub_assembly)
-        info("Assembly: ", round(j/nchunks*100,1), " % done. ")
+        if ne > 100
+            info("Assembly: ", round(j/nchunks*100,1), " % done. ")
+        end
     end
 #   optimize!(assembly)
 #   dim = length(assembly.stiffness_matrix.I)

src/elements.jl

@@ -210,3 +210,26 @@ function Base.haskey(element::Element, what)
     haskey(element.fields, what)
 end
 
+""" Calculate local normal-tangential coordinates for element. """
+function calculate_normal_tangential_coordinates!{E}(element::Element{E}, time::Real)
+    proj(u, v) = dot(v, u) / dot(u, u) * u
+    ntcoords = Matrix[]
+    refcoords = get_reference_element_coordinates(E)
+    x = element("geometry", time)
+    for xi in refcoords
+        dN = get_dbasis(E, xi)*x
+        normal = cross(dN[:,1], dN[:,2])
+        normal /= norm(normal)
+        u1 = normal
+        j = indmax(abs(u1))
+        v2 = zeros(3)
+        v2[mod(j,3)+1] = 1.0
+        u2 = v2 - proj(u1, v2)
+        u3 = cross(u1, u2)
+        tangent1 = u2/norm(u2)
+        tangent2 = u3/norm(u3)
+        push!(ntcoords, [normal tangent1 tangent2])
+    end
+    element["normal-tangential coordinates"] = ntcoords
+end
+

src/lagrange.jl

@@ -27,6 +27,10 @@ function calculate_lagrange_basis_coefficients(P, X)
     return inv(A)'
 end
 
+function refcoords(X::Matrix)
+    return Vector{Float64}[X[:,i] for i=1:size(X,2)]
+end
+
 """
 Create new Lagrange element
 
@@ -37,7 +41,10 @@ Examples
 macro create_lagrange_element(element_name, element_description, X, P)
     eltype = esc(element_name)
     quote
-        global get_basis, get_dbasis
+        global get_basis, get_dbasis,
+               get_reference_element_coordinates,
+               get_reference_element_midpoint
+
         #basis, dbasis = calculate_lagrange_basis($P, $X)
         C = calculate_lagrange_basis_coefficients($P, $X)
         basis(xi) = C*$P(xi)
@@ -49,6 +56,15 @@ macro create_lagrange_element(element_name, element_description, X, P)
             return basis(xi)'
         end
 
+        XX = refcoords($X)
+        function get_reference_element_coordinates(::Type{$eltype})
+            return XX
+        end
+
+        XXX = vec(mean($X, 2))
+        function get_reference_element_midpoint(::Type{$eltype})
+            return XXX
+        end
 #=
         function get_dbasis(::Type{$eltype}, xi::Vector{Float64})
             return dbasis(xi)'
@@ -81,6 +97,10 @@ end
      0.0 0.0 1.0],
     (xi) -> [1.0, xi[1], xi[2]])
 
+#function get_reference_element_midpoint(::Type{Tri3})
+#    return [1.0/3.0, 1.0/3.0]
+#end
+
 @create_lagrange_element(Tri6, "6 node quadratic triangle element",
     [0.0 1.0 0.0 0.5 0.5 0.0
      0.0 0.0 1.0 0.0 0.5 0.5],

src/mortar.jl

@@ -11,7 +11,7 @@ function project_from_slave_to_master{S,M}(slave::Element{S}, master::Element{M}
 
     # slave side geometry and normal direction at xi1
     X1 = slave("geometry", xi1, time)
-    N1 = slave("nodal ntsys", xi1, time)[:,1]
+    N1 = slave("normal-tangential coordinates", xi1, time)[:,1]
 
     # master side geometry at xi2
     #master_basis = master.basis.data.basis
@@ -59,7 +59,7 @@ function project_from_master_to_slave{S,M}(slave::Element{S}, master::Element{M}
     # slave side geometry and normal direction at xi1
 
     slave_geometry = slave("geometry")(time)
-    slave_normals = slave("nodal ntsys")(time)
+    slave_normals = slave("normal-tangential coordinates")(time)
     #slave_basis = slave.basis.data.basis
     #slave_dbasis = slave.basis.data.dbasis
     slave_basis(xi) = get_basis(S, [xi])
@@ -86,7 +86,7 @@ function project_from_master_to_slave{S,M}(slave::Element{S}, master::Element{M}
     end
 
     #X1(xi1) = slave_basis("geometry", [xi1], time)
-    #N1(xi1) = slave_basis("nodal ntsys", [xi1], time)[:,1]
+    #N1(xi1) = slave_basis("normal-tangential coordinates", [xi1], time)[:,1]
 
     #master_basis = get_basis(master)
 
@@ -187,11 +187,16 @@ Notes
 
 """
 # function create_auxiliary_plane(x, ximp, normals, basis)
-function create_auxiliary_plane(element::Element{Tri3}, time::Real)
-    proj(u, v) = dot(v, u) / dot(u, u) * u
-    xi = [1.0/3.0, 1.0/3.0]
+function create_auxiliary_plane{E}(element::Element{E}, time::Real)
+#   proj(u, v) = dot(v, u) / dot(u, u) * u
+#    xi = [1.0/3.0, 1.0/3.0]
+
+    xi = get_reference_element_midpoint(E)
     x0 = element("geometry", xi, time)
-    n = element("nodal ntsys", xi, time)[:, 1]
+    ntbasis = element("normal-tangential coordinates", xi, time)
+    return x0, ntbasis
+#=
+    n = element("normal-tangential coordinates", xi, time)[:, 1]
     n /= norm(n)
     # gram-schmidt
     u1 = n
@@ -204,6 +209,7 @@ function create_auxiliary_plane(element::Element{Tri3}, time::Real)
     t2 = u3/norm(u3)
     new_basis = [n t1 t2]
     return x0, new_basis
+=#
 end
 
 """
@@ -243,7 +249,14 @@ function project_point_to_auxiliary_plane(p::Vector, x0::Vector, Q::Matrix)
     n = Q[:,1]
     ph = p - dot(p-x0, n)*n
     qproj = Q'*(ph-x0)
-    @assert isapprox(qproj[1], 0.0)
+    if !isapprox(qproj[1], 0.0; atol=1.0e-12)
+        info("project_point_to_auxiliary_plane(): point not projected correctly.")
+        info("p: $p")
+        info("x0: $x0")
+        info("Q: \n$Q")
+        info("qproj: $qproj")
+        error("Failed to project point to auxiliary plane.")
+    end
     return qproj[2:3]
 end
 
@@ -626,6 +639,7 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly, problem::BoundaryPro
             N2 = master_element(xi_projected, time)
             S = w*N1'*N1
             M = w*(N1'*N2)'
+#           M = w*N1'*N2
             # FIXME: why this needs now to be transpose?
             # assembly / repeat
             for i=1:field_dim