Updated concept of elements.

src/JuliaFEM.jl

@@ -2,9 +2,10 @@
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 module JuliaFEM
 
-
 VERSION < v"0.4-" && using Docile
 using Lexicon
+using Logging
+@Logging.configure(level=DEBUG)
 
 include("types.jl") # type definitions
 include("math.jl") # basic mathematical operations

src/abaqus_reader.jl

@@ -1,18 +1,11 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-module abaqus_reader
-
-using Logging
-@Logging.configure(level=DEBUG)
-
-VERSION < v"0.4-" && using Docile
-
 eldims = Dict(
     "C3D10" => 10,
     "C3D4" => 4)
-global handlers = Dict()
 
+global handlers = Dict()
 
 """
 Register new handler for parser
@@ -57,7 +50,7 @@ function parse_element_section(model, header, data)
     end
     eldim = eldims[eltype]
     m = matchall(r"[0-9]+", data)
-    m = map(integer, m)
+    m = map((s) -> parse(Int, s), m)
     elements = create_or_get(model, "elements")
     m = reshape(m, eldim+1, round(Int, length(m)/(eldim+1)))
     nel = size(m)[2]
@@ -81,7 +74,7 @@ function parse_nodeset_section(model, header, data)
     nset_name = header["options"]["NSET"]
     Logging.debug("Creating node set $nset_name")
     m = matchall(r"[0-9]+", data)
-    node_ids = map(integer, m)
+    node_ids = map((s) -> parse(Int, s), m)
     nsets = create_or_get(model, "nsets")
     nsets[nset_name] = Int64[]
     for j in node_ids
@@ -132,5 +125,3 @@ add_handler("NODE", parse_node_section)
 add_handler("ELEMENT", parse_element_section)
 add_handler("NSET", parse_nodeset_section)
 
-end
-

src/math.jl

@@ -5,9 +5,10 @@
 This module contains math stuff, including interpolation, integration, linearization, ...
 """
 
-using JuliaFEM
 using ForwardDiff
 
+export interpolate, integrate, linearize
+
 """
 Interpolate field variable using basis functions f for point ip.
 This function tries to be as general as possible and allows interpolating
@@ -59,7 +60,7 @@ function interpolate{T<:Real}(field::Array{T,2}, basis::Function, ip)
     return result
 end
 function interpolate(e::Element, field::ASCIIString, x::Array{Float64,1}; derivative=false)
-    return interpolate(e.attributes[field], derivative ? e.dbasis : e.basis, x)
+    return interpolate(e.attributes[field], derivative ? e.shape_functions.dbasis : e.shape_functions.basis, x)
 end
 
 
@@ -67,15 +68,15 @@ end
 
 """
 function get_basis(el::Element, xi)
-    return el.basis(xi)
+    return el.shape_functions.basis(xi)
 end
 
 """
 Return partial derivatives of shape functions w.r.t X using chain rule.
 """
-function get_dbasisdX(el::Element, ip)
-    J = interpolate(el, "coordinates", ip.xi; derivative=true)
-    dbasisdX = el.dbasis(ip.xi)*inv(J')
+function get_dbasisdX(el::Element, xi)
+    J = interpolate(el, "coordinates", xi; derivative=true)
+    dbasisdX = el.shape_functions.dbasis(xi)*inv(J')
     return dbasisdX
 end
 
@@ -96,7 +97,7 @@ Array{Float64, 2}
     jacobian / "tangent stiffness matrix"
 
 """
-function linearize(f::Function, el::JuliaFEM.Element, field::ASCIIString)
+function linearize(f::Function, el::Element, field::ASCIIString)
     dim, nnodes = size(el.attributes[field])
     function helper!(x, y)
         orig = copy(el.attributes[field])
@@ -112,7 +113,7 @@ end
 This version returns another function which can be then evaluated against field
 """
 function linearize(f::Function, field::ASCIIString)
-    function jacobian(el::JuliaFEM.Element, args...)
+    function jacobian(el::Element, args...)
         dim, nnodes = size(el.attributes[field])
         function helper!(x, y)
             orig = copy(el.attributes[field])
@@ -129,7 +130,7 @@ end
 """
 In-place version, no additional garbage collection.
 """
-function linearize!(f::Function, el::JuliaFEM.Element, field::ASCIIString, target::ASCIIString)
+function linearize!(f::Function, el::Element, field::ASCIIString, target::ASCIIString)
     el.attributes[target][:] = 0.0
     dim, nnodes = size(el.attributes[field])
     function helper!(x, y)
@@ -154,23 +155,31 @@ el::Element
 f::Function
     Function to integrate
 """
-function integrate(f::Function, el::JuliaFEM.Element)
+function integrate(f::Function, el::Element)
     target = []
     for ip in el.integration_points
-        J = JuliaFEM.interpolate(el, "coordinates", ip.xi; derivative=true)
+        J = interpolate(el, "coordinates", ip.xi; derivative=true)
         push!(target, ip.weight*f(el, ip)*det(J))
     end
     return sum(target)
 end
+#function integrate(f::Function, integration_points::Array{IntegrationPoint, 1}, Xargs...)
+#    target = []
+#    for ip in integration_points
+#        J = interpolate(el, "coordinates", ip.xi; derivative=true)
+#        push!(target, ip.weight*f(ip, args...)*det(J))
+#    end
+#    return sum(target)
+#end
 
 """
 This version returns a function which must be operated with element e
 """
 function integrate(f::Function)
-    function integrate(el::JuliaFEM.Element)
+    function integrate(el::Element)
         target = []
         for ip in el.integration_points
-            J = JuliaFEM.interpolate(el, "coordinates", ip.xi; derivative=true)
+            J = interpolate(el, "coordinates", ip.xi; derivative=true)
             push!(target, ip.weight*f(el, ip)*det(J))
         end
         return sum(target)
@@ -181,11 +190,12 @@ end
 """
 This version saves results inplace to target, garbage collection free
 """
-function integrate!(f::Function, el::JuliaFEM.Element, target)
+function integrate!(f::Function, el::Element, target)
     # set target to zero
     el.attributes[target][:] = 0.0
     for ip in el.integration_points
-        J = JuliaFEM.interpolate(el, "coordinates", ip.xi; derivative=true)
+        J = interpolate(el, "coordinates", ip.xi; derivative=true)
         el.attributes[target][:,:] += ip.weight*f(el, ip)*det(J)
     end
 end
+

src/types.jl

@@ -1,6 +1,8 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
+export IntegrationPoint, Element, Assembly, FunctionSpace
+
 """
 Integration point
 
@@ -18,18 +20,21 @@ type IntegrationPoint
     attributes :: Dict{ASCIIString, Any}
 end
 
-type Element
-    id :: Int
-#   element_type :: Int
-    node_ids :: Array{Int, 1}
+type FunctionSpace
     basis :: Function
     dbasis :: Function
-    integration_points :: Array{IntegrationPoint, 1}
-    attributes :: Dict{ASCIIString, Any}
-#   ipoints :: Array{Float64, 2}
-#   iweights :: Array{Float64, 1}
 end
 
+abstract Element
+
+#type Element
+#    id :: Int
+#    node_ids :: Array{Int, 1}
+#    shape_functions :: FunctionSpace
+#    integration_points :: Array{IntegrationPoint, 1}
+#    attributes :: Dict{ASCIIString, Any}
+#end
+
 
 type Assembly
     # LHS
@@ -42,4 +47,6 @@ type Assembly
     # global dofs for each element
     gdofs :: Dict{Int64, Array{Int64, 1}}
 end
+Assembly() = Assembly(Int64[], Int64[], Float64[], Int64[], Float64[], Dict{Int64,Array{Int64,1}}())
+Assembly(gdofs::Dict{Int64,Array{Int64,1}}) = Assembly(Int64[], Int64[], Float64[], Int64[], Float64[], gdofs)
 

src/xdmf.jl

@@ -1,17 +1,8 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-module xdmf
-
-using Logging
-@Logging.configure(level=INFO)
-
 using LightXML
 
-VERSION < v"0.4-" && using Docile
-
-# i add docstrings later
-
 # element codes: http://www.paraview.org/pipermail/paraview/2013-July/028859.html
 # > from  ./VTK/ThirdParty/xdmf2/vtkxdmf2/libsrc/XdmfTopology.h
 # >
@@ -75,29 +66,68 @@ function xdmf_new_grid(temporal_collection; time=0)
     return grid
 end
 
+#function xdmf_new_mesh(grid, X, elmap)
+#    geometry = new_child(grid, "Geometry")
+#    set_attribute(geometry, "Type", "XYZ")
+#    dataitem = new_child(geometry, "DataItem")
+#    set_attribute(dataitem, "DataType", "Float")
+#    set_attribute(dataitem, "Dimensions", length(X))
+#    set_attribute(dataitem, "Format", "XML")
+#    set_attribute(dataitem, "Precision", "4")
+#    add_text(dataitem, join(X, " "))
+#    topology = new_child(grid, "Topology")
+#    set_attribute(topology, "Dimensions", "1")
+#    set_attribute(topology, "Type", "Mixed")
+#    dataitem = new_child(topology, "DataItem")
+#    set_attribute(dataitem, "DataType", "Int")
+#    set_attribute(dataitem, "Dimensions", length(elmap))
+#    set_attribute(dataitem, "Format", "XML")
+#    set_attribute(dataitem, "Precision", 4)
+#    elmap2 = copy(elmap)
+#    elmap2[2:end,:] -= 1
+#    add_text(dataitem, join(elmap2, " "))
+#end
+
 function xdmf_new_mesh(grid, X, elmap)
+    dim, nnodes = size(X)
     geometry = new_child(grid, "Geometry")
     set_attribute(geometry, "Type", "XYZ")
     dataitem = new_child(geometry, "DataItem")
     set_attribute(dataitem, "DataType", "Float")
-    set_attribute(dataitem, "Dimensions", length(X))
+    set_attribute(dataitem, "Dimensions", "$nnodes $dim")
     set_attribute(dataitem, "Format", "XML")
-    set_attribute(dataitem, "Precision", "4")
-    add_text(dataitem, join(X, " "))
+    set_attribute(dataitem, "Precision", 8)
+    #add_text(dataitem, join(X, " "))
+    s = "\n"
+    
+    for i=1:nnodes
+        s *= "\t\t" * join(X[:,i], " ") * "\n"
+    end
+    s *= "       "
+    add_text(dataitem, s)
 
-    topology = new_child(grid, "Topology")
-    set_attribute(topology, "Dimensions", "1")
-    set_attribute(topology, "Type", "Mixed")
-    dataitem = new_child(topology, "DataItem")
-    set_attribute(dataitem, "DataType", "Int")
-    set_attribute(dataitem, "Dimensions", length(elmap))
-    set_attribute(dataitem, "Format", "XML")
-    set_attribute(dataitem, "Precision", 4)
     elmap2 = copy(elmap)
     elmap2[2:end,:] -= 1
-    add_text(dataitem, join(elmap2, " "))
+    dim, nelements = size(elmap2)
+
+    topology = new_child(grid, "Topology")
+    #set_attribute(topology, "Dimensions", "1")
+    set_attribute(topology, "TopologyType", "Mixed")
+    set_attribute(topology, "NumberOfElements", nelements)
+    dataitem = new_child(topology, "DataItem")
+    set_attribute(dataitem, "DataType", "Int")
+    set_attribute(dataitem, "Dimensions", "$nelements $dim")
+    set_attribute(dataitem, "Format", "XML")
+    set_attribute(dataitem, "Precision", 8)
+    s = "\n"
+    for i=1:nelements
+        s *= "\t\t" * join(elmap2[:,i], " ") * "\n"
+    end
+    add_text(dataitem, s)
+    #add_text(dataitem, join(elmap2, " "))    
 end
 
+
 function xdmf_new_field(grid, name, source, data)
     loc = Dict("elements" => "Cell",
                "nodes" => "Node")
@@ -140,4 +170,3 @@ function xdmf_save_model(xdoc, filename)
   save_file(xdoc, filename)
 end
 
-end

test/test_abaqus_reader.jl

@@ -5,45 +5,61 @@ using FactCheck
 using Logging
 @Logging.configure(level=INFO)
 
-using JuliaFEM.abaqus_reader: parse_abaqus, parse_element_section
+#using JuliaFEM.abaqus_reader: parse_abaqus, parse_element_section
+include(Pkg.dir("JuliaFEM")*"/src/abaqus_reader.jl")
 
 facts("test import abaqus model") do
   # FIXME: get_test_data()
     fid = open(Pkg.dir("JuliaFEM")*"/geometry/3d_beam/palkki.inp")
     model = parse_abaqus(fid)
     close(fid)
-    @fact length(model["nodes"]) => 298
-    @fact length(model["elements"]) => 120
-    @fact length(model["elsets"]["Body1"]) => 120
-    @fact length(model["nsets"]["SUPPORT"]) => 9
-    @fact length(model["nsets"]["LOAD"]) => 9
-    @fact length(model["nsets"]["TOP"]) => 83
+    @fact length(model["nodes"]) --> 298
+    @fact length(model["elements"]) --> 120
+    @fact length(model["elsets"]["Body1"]) --> 120
+    @fact length(model["nsets"]["SUPPORT"]) --> 9
+    @fact length(model["nsets"]["LOAD"]) --> 9
+    @fact length(model["nsets"]["TOP"]) --> 83
 end
 
 facts("test that reader throws error when dimension information of elemenet is missing") do
-#   *ELEMENT, TYPE=neverseenbefore, ELSET=Body1
-  data = """
-         1,       243,       240,       191,       117,       245,       242,       244,
-         1,         2,       196
-  """
-  model = Dict()
-  header = Dict("section"=>"ELEMENT", "options" => Dict("TYPE" => "neverseenbefore", "ELSET"=>"Body1"))
-  @fact_throws parse_element_section(model, header, data)
+    #   *ELEMENT, TYPE=neverseenbefore, ELSET=Body1
+    data = """
+    1,       243,       240,       191,       117,       245,       242,       244,
+    1,         2,       196
+    """
+    model = Dict()
+    header = Dict("section"=>"ELEMENT", "options" => Dict("TYPE" => "neverseenbefore", "ELSET"=>"Body1"))
+    @fact_throws parse_element_section(model, header, data)
+end
+
+facts("read element section") do
+    data = """
+    1,       243,       240,       191,       117,       245,       242,       244,
+    1,         2,       196
+    2,       204,       199,       175,       130,       207,       208,       209,
+    3,         4,       176
+    """
+    model = Dict()
+    header = Dict("section" => "ELEMENT", "options" => Dict("TYPE" => "C3D10", "ELSET" => "BEAM"))
+    parse_element_section(model, header, data)
+    @fact length(model["elements"]) --> 2
+    @fact model["elements"][1] --> [243, 240, 191, 117, 245, 242, 244, 1, 2, 196]
+    @fact model["elements"][2] --> [204, 199, 175, 130, 207, 208, 209, 3, 4, 176]
 end
 
 facts("test unknown handler warning message") do
-  fn = tempname()
-  fid = open(fn, "w")
-  testdata = """
-  *ELEMENT2, TYPE=C3D10, ELSET=Body1
-         1,       243,       240,       191,       117,       245,       242,       244,
-         1,         2,       196
-  """
-  write(fid, testdata)
-  close(fid)
-  fid = open(fn)
-  model = parse_abaqus(fid)
-  close(fid)
-  # empty model expected, parser doesn't know what to do with unknown section
-  @fact length(model) => 0
+    fn = tempname()
+    fid = open(fn, "w")
+    testdata = """
+    *ELEMENT2, TYPE=C3D10, ELSET=Body1
+    1,       243,       240,       191,       117,       245,       242,       244,
+    1,         2,       196
+    """
+    write(fid, testdata)
+    close(fid)
+    fid = open(fn)
+    model = parse_abaqus(fid)
+    close(fid)
+    # empty model expected, parser doesn't know what to do with unknown section
+    @fact length(model) --> 0
 end